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THE  NEC2  RADIATION  PATTERNS  OF  UNDER-SEGMENTED  WIRE  GRID  MODELS  OF 
A  FIGHTER  AIRCRAFT  COMPARED  TO  MEASUREMENTS 

O.  Givati  and  APC.  Fourie 
University  of  the  Witwatersrand,  Johannesburg 
Private  Bag  3,  P  O  WITS  2050,  South  Africa 


ABSTRACT 

Hiis  paper  presents  results  of  a  NEC2  method  of  moments 
evaluation  of  VHF/UHF  antennas  on  a  wire  grid  model  of 
a  frghter  aircraft  The  study  shows  that  practically  usefril 
radiation  patterns  results  can  be  obtained  when  the  grid 
model  is  considerably  under-segmented  in  aircraft  regions 
which  are  electrically  far  removed  from  the  antennas. 
Normal  modelling  guide-lines  requites  approximately  35000 
segments  for  the  fighter  grid  model  at  400  MHz;  judicious 
under-segmentation  gave  useful  results  using  a  model 
comprising  5000  segments.  This  reduction  in  segments 
reduced  computer  time  by  a  factor  of  343  and  memory 
requirements  is  reduced  by  a  factor  of  18.  NEC2  radiation 
patterns  are  compared  to  measurements  on  a  1/lOth  scale 
model  which  was  performed  in  an  anechoic  chamber 
compact  range. 

1  INTRODUCTION 

The  paper  consider  a  stucfy  of  radiation  patterns  from 
antennas  mounted  on  a  fighter  aircraft  NEC2  [1]  computed 
radiation  patterns  are  used  for  subsequent  statistical 
assessment  of  link  performance  for  the  aircraft  during  typical 
mission  profiles.  Wire  grids  and  segments  much  longer  than 
the  recommended  [1]  O.IA,  were  used  in  areas  far  removed 
from  the  anteimas  to  reduce  execution  time  and  memory 
requirements.  Measured  results,  as  well  as  comparison  to 
more  densely  segmented  numerical  modeb,  were  used  to 
ensure  that  these  violations  of  modelling  rules  still  produced 
usefril  values.  It  may  be  argued  that  numerical  modelling  is 
superfluous  when  measured  resulte  are  available,  but  this  is 
not  the  case  when  considering  that: 

*  measuremenb  in  only  three  principle  planes  wete 
performed  and  compared  to  a  subset  of  the  computed 
values;  the  statbtical  link  analysis  requites  full  three 
dimensional  pattern  information  which  is  difficult  and 
time  consuming  to  obtain  by  measurement 

•  traditional  engineering  normally  argues  that 
measuremenb  constitute  the  more  defiiutive 
characterization  of  a  system  when  compared  with 
calculation.  This  is  definitely  not  always  the  case  with 


radiation  pattern  measuremenb,  especially  on  scale 
modeb,  and  many  electromagnetics  engineers  can  attest 
to  cases  where  more  confidence  can  be  pbced  on  values 
obbined  from  calculation  or  simulation. 

The  engineer  b  hence  in  a  difficult  position.  He  has  two 
methods  giving  two  seb  of  resulb  -  both  associated  with 
some  potentially  large  errors.  On  the  positive  side  b  the  fact 
that  the  two  seb  of  resulb  were  obbined  using  two  entirely 
difrerent  methods,  with  appropriate  techruques  used  in  both 
cases  to  minimize  errors.  Qualibtive  agreement  will 
definitely  demonstrate  the  absence  of  major  blunders  in  both 
methods,  and  engineering  judgement  and  cost  will  dicbte 
which  set  of  resulb  is  likely  to  be  more  appropriate. 

Valid  numerical  modeb,  for  the  purpose  of  subsequent 
sbtbtical  link  analysis,  are  hence  those  which  show 
qualibtive  agreement  with  measured  resulb.  Quantibtive 
comparison  between  the  measured  and  computed  radiation 
patterns  is  difficult  and  often  misleading,  because: 

•  errors  may  be  large  when  comparing  measured  and 
computed  values  at  specific  angles,  but  such  errors  may 
only  be  due  to  a  slight  offset  in  the  position  of  a 
radiation  pattern  null,  for  instance.  Such  small  angular 
offseb  are  not  of  any  concern  when  performing 
sbtbtical  analysis  of  communication  link  performance. 

*  the  existence  of  significant  measurement  errors  also 
frustrate  efforb  to  call  the  deviation  from  measured 
values  "simulation  errors". 

It  is  hence  evident  that  some  of  the  more  traditional 
quantibtive  measures  of  agreement  between  measured  and 
computed  resulb  will  be  less  useful  and  often  meaningless. 
It  is  indeed  left  to  the  reader  to  assess  the  presented 
comparisons  between  measured  and  computed  resulb  and 
decide  on  their  worth  for  a  specific  application.  The  most 
difficult  aspect  of  the  comparison,  in  fact,  is  the  inaccuracies 
with  the  measuremenb  themselves;  the  extent  of  such 
inaccuracies  can  easily  be  gouged  by  the  deviations  from 
symmetry  in  certain  planes,  as  well  as  comparing 
corresponding  poinb  where  different  planes  intersect 
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2  SCALE  MODEL  MEASUREMENTS 

A  1/1 0th  scale  model  of  the  fighter  aircraft  was  constructed 
from  polystyrene  and  covered  with  polymer  sheets.  The 
model  was  then  copper-plated  to  ensure  adequate 
conductivity.  The  antennas  were  incorporated  into  the  model 
by  intemaUy  guiding  coaxial  cable  to  the  antenna  positions 
and  protruding  the  inner  conductor  by  30mm  to  foim 
monopoles  corresponding  to  the  top-fin  and  bottom-fin 
antennas  respectively. 

Radiation  pattern  measurements  in  the  principle  planes 
(azimuth,  side  and  pitch  roll)  were  then  performed  on  the 
scale  model  using  a  compact  range  suitable  for  the  frequency 
range  2GHz-18GHz.  The  compact  range  consists  of  a 
25  m  X  10  m  X  10  m  anechoic  chamber,  calibrated  feed 
antenna,  offset  parabolic  reflector  and  positioner  with  three 
axes  of  freedom.  The  range  was  calibrated  using  a  vertically 
polarized,  standard  gain  horn  in  the  frequency  range 
2  GHz-4  GHz.  The  standard  gain  antenna  was  then  replaced 
with  the  scale  model  with  the  mounting  bar  in  the  .same 
horizontal  position  as  the  standard  gain  antenna  and  the 
middle  of  the  model  at  the  same  height  as  the  standard  gain 
antenna.  The  measurement  system  hence  measured  absolute 
gain  in  dB  relative  to  an  isotropic  source  (dBi).  The  only 
adjustment  to  the  measured  dBi  values  was  to  compute  the 
losses  in  the  cable  leading  from  the  calibrated  connector  to 
the  antenna  and  to  subtract  these  losses  from  the 
measurements  in  order  to  obtain  the  actual  gain  as  measured 
at  the  antenna  port 

The  following  movements  were  executed  with  the  positioner 
to  measure  the  pattern  in  the  three  principle  aircraft  planes 
(only  the  position  and  movements  for  the  top  fin  antenna  are 
given  below;  the  bottom  fin  was  measured  using  exact 
inverse  positions  and  positioner  movements): 

•  Azimuth  (yaw)  plane:  The  model  was  mounted 
horizontally  and  upright  and  simply  rotated  through 
360°. 

•  Pitch  plane:  The  model  was  mounted  upright  and 
horizontal,  facing  the  receiving  antenna.  The  model  was 
then  tilted  90°  forward  and  90°  backwards  to  produce 
the  measurement  points  between  0°  and  180°  in  the  plots 
of  results.  The  model  was  then  mounted  facing  away 
from  the  receiving  antenna  and  again  tilted  90°  towards 
and  90°  away  from  the  receiving  antenna  to  produce  the 
values  plotted  between  180°  through  360°.  It  is  during 
the  extremes  of  these  tilting  movements  away  from  the 
receiving  antenna  that  reflections  associated  with  the 


positioner  and  mounting  bar  bending  (referred  to  below) 
were  most  severe  (corresponding  to  the  plotted  values 
at  angles  0°  and  180°). 

•  Roll  plane:  The  model  was  mounted  upright  and 
broadside  to  the  receiving  antenna  (one  wing  pointing 
towards  the  receiving  antenna).  The  model  was  then 
tilted  90°  towards  and  90°  away  from  the  receiving 
antenna  to  produce  one  half  of  the  measured  points.  It 
was  then  rotated  by  180°  and  the  same  movements  were 
repeated  to  yield  the  remaining  values  of  measurements. 
Once  again  the  comments  above  regarding  errors  apply 
with  most  severe  cases  for  roll  plane  graph  angels  90° 
and  270°. 

It  should  be  noted  that  only  vertically  polarized  gain  was 
measured,  since  this  is  the  dominant  polarization  from  both 
antennas  and  was  also  the  only  polarization  of  interest  when 
performing  link  assessments.  The  signal  source  for 
measurements  was  also  linked  to  the  scale  model  via  a  cable 
which  was  always  routed  along  the  body  of  the  scaled  model 
to  ensure  that  the  cable  enters  the  model  at  the  opposite  side 
of  the  model  in  relation  to  the  antenna  position  (model 
mounting  upright  and  upside  down  was  possible  using  the 
mounting  arrangement  to  facilitate  this  aim).  When  HF 
measurements  on  scale  models  were  performed  in  the  past, 
a  stand  alone  source  was  constructed  and  housed  inside  the 
scale  model  fuselage,  because  at  lower  frequencies  the 
aircraft  electrical  dimensions  are  small  in  terms  of 
wavelengths,  and  the  cable  shield  carries  substantial  currents 
which  affect  the  measurements.  At  higher  frequencies, 
however,  the  interaction  between  an  antenna  mounted  on  one 
side  of  the  aircraft  with  the  measurement  cable  on  the 
opposite  side  is  minimal.  A  stand  alone  signal  source  would 
have  been  a  disadvantage,  in  the  VHF/UHF  case  where  one 
is  interested  in  absolute  gain  values,  because  it  would  need 
to  be  custom  designed  with  suitable  calibrated 
characteristics. 

No  measurements  are  error-free,  and  measurement 
uncertainty  is  particularly  difficult  to  ascertain  when 
radiation  patterns  are  measured.  Using  the  described 
measurement  set-up  the  following  factors  may  have  been  the 
cause  of  some  of  the  errors  (these  are  more  or  less  listed  in 
order  of  their  severity): 

•  Inaccuracies  in  the  scale  model  finish  and  dimensions. 
Some  of  the  results  presented  shows  some  signs  of 
asymmetry  which  is  most  likely  due  to  slight  errors  in 
curvature  on  either  side  of  the  antenna 
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■  Literacdon  between  the  model  and  the  positionet.  These 
ate  most  severe  when  the  model  is  tilted  backwards 
(away  from  the  source)  during  roll  and  pitch  plane 
measurements.  This  interaction  manifests  itself  most 
prominently  when  one  observes  the  errors  in  the 
measurements  at  the  angles  where  the  measured  values 
do  not  coincide.  (Note  the  90°  and  270°  measurements 
in  roll  plane  and  the  0°  and  180°  points  during  pitch 
plane  measurements). 

•  Angular  errors  due  to  a  certain  amount  of  bending  in 
the  perspex  mounting  bar  when  the  model  was  tilted  for 
roll  and  pitch  plane  measurements.  Once  again  most 
severe  for  the  90°  and  270°  measurements  in  roll  plane 
and  the  0°  and  180°  points  during  pitch  plane 
measurements. 

•  Rotation  "wobbles"  because  the  mount  in  the  fuselage 
is  not  absolutely  perpendicular  to  the  fuselage  horizontal 
datum  line. 

•  Errors  due  to  the  fact  that  the  protruding  monopoles 
were  not  always  exactly  straight.  This  manifests  itself 
most  commonly  in  terms  of  errors  in  the  position  of  the 
natural  monopole  pattern  nulls. 

No  significant  errors  were  caused  by  source  amplitude 
instability,  reflecdons  from  the  chamber  walls  and  other 
^stematic  errors  associated  with  the  range.  The  compact 
range  used  is  professionally  constructed  and  designed  for 
measurements  in  the  frequency  range  of  interest.  The 
calibration  procedure  ensures  that  equipment  error  levels  are 
accounted  for,  and  the  quiet  zone  associated  with  the 
measurements  is  larger  than  the  maximum  movement  and 
model  dimensions  during  the  course  of  the  measurement 
manoeuvres. 

3  GENERATION  OF  THE  NUMERICAL 
MODELS 

In  this  section,  representations  of  the  aircraft  grid  models, 
used  for  the  evaluation  of  the  top  and  bottom  fin  antennas’ 
performance,  are  shown  in  Fig.  4  and  6.  Different  grid 
models  were  constructed  in  order  to  evaluate  the  antermas’ 
performances  at  different  frequencies.  The  numerical  model 
uses  symmetry,  and  therefore  only  half  of  each  representative 
grid  model  is  displayed.  The  two  halves  are  coimected  at  the 
aircraft  axis  of  symmetry,  where  the  segments  are  seen  to 
be  terminated. 


The  numerical  models  of  the  aircraft  were  constructed  using 
the  Structure  Interpolation  and  Gridding  software  package, 
SIG  [2],  developed  by  EM-Simulations  (Pty)  Ltd.  The  SIG 
package  generates  a  three-dimensional  grid  model  from  a  set 
of  user  defined  cross-sectional  cuts  at  points  of  abrupt 
change  along  the  three  dimensional  structure  -  as  indicated 
in  Fig.  3.  The  user-defined  cross-sections  shown  in  Fig.  3 
form  the  basis  for  the  grid  models,  shown  in  Figs.  4  and  S. 
Wings  and  other  features  attached  to  the  main  fuselage  are 
accommodated,  using  the  SIG  package,  by  tagging  curves 
which  represents  features  in  a  user  defined  cross  section  and 
using  corresponding  tag  numbers  in  later  user  defined  cross 
sections.  In  this  way,  for  instance,  cross-section  11  in  Fig. 
3  consists  of  3  curves  with  three  tag  numbers:  The  first  curve 
represents  the  main  fuselage,  the  second  the  top  extension 
of  the  cockpit  and  the  last  curve  in  that  cross  section  will 
just  be  a  single  point  representing  the  start  of  the  wing.  In 
the  next  cross  section  (12  in  Fig.  3)  these  three  curves  are 
re-defined  with  the  dimensions  at  that  cross  sectional  point 
and  the  curve  which  was  tagged  as  the  wing  will  in  this  case 
be  a  line;  interpolation  between  the  single  point  curve  in 
cross  section  11  and  the  line  in  cross  section  12  defines  the 
characteristic  delta  wing  of  this  fighter  aircraft  The  ability 
of  SIG  to  accommodate  appendages  to  a  fuselage  in  this 
fashion  is  exceedingly  useful,  because  the  gridding  routine 
ensures  that  the  attached  features  are  connected  at  all  points 
of  the  grid. 

The  grid  models  shown  in  Figs.  4  and  5  are  generated  by 
interpolating  the  cross-sectional  cuts  between  the  user 
defined  cross-sections  (in  Fig.  3),  at  intervals  which  are  not 
greater  than  the  specified  target  segment  lengths  of  the 
user-defined  cross-sections.  The  segmentations  produced  by 
SIG  are  mainly  quadrilaterals,  with  the  side  lengths 
approximately  equal  to  the  target  segment  lengths  requested. 
Some  triangular  grid  elements  are  also  formed  when  curves 
expand  or  contract  from  one  cross-section  to  the  next  (eg. 
wings  regions).  These  triangular  grid  elements  also  have 
edge  lengths  approximately  equal  to  the  required  segment 
length.  The  segment  radii  are  calculated  by  the  SIG  package 
to  ensure  that  the  surface  area  of  the  segm^ts  comprising 
the  grid  is  approximately  twice  the  surface  area  of  the 
structure  which  is  modelled.  Also  note  that  SIG 
automatically  generated  segments  abutting  the  ^mmetry 
plane  which  are  only  half  the  grid  length  to  ensure  grid  size 
continuity  across  the  synunetiy  plane. 
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The  SIG  program  allows  the  user  to  specify  a  specific  target 
segment  length  for  every  user  deflned  cross  section.  This 
target  segment  length  is  then  used  until  a  new  target  length 
is  specified  in  a  subsequent  user  defined  cross  section.  This 
SIG  feature  allowed  the  variable  segmentation  (and  hence 
grid  size)  used  during  this  study. 

The  use  of  segment  lengths  larger  than  0.1  A,  significantly 
reduces  the  number  of  segments  used  in  the  numerical 
model.  Generally,  such  a  violation  of  the  numerical 
modelling  rules  can  result  in  an  invalid  numerical  model.  In 
this  paper,  we  show  that  the  effect  of  imder-segmentation, 
in  regions  which  are  sufficiently  removed  from  the  antermas, 
is  negligible  on  electrically  large  structures.  In  return,  the 
reduction  in  computation  time  is  significantly  large.  The 
effect  of  under-segmentation  can  be  seen  in  Fig.  1  and  2, 
where  the  radiation  patterns  of  the  top-fin  antetma  on  two 
different  120MHz  grid  models  is  simulated  according  to  the 
segmentation  schemes  indicated  in  Table  1. 


Section 

No. 

(Fig-  3) 

Target 

segment 

length 

(120MH2) 

Target 

segment 

length 

(120MHz) 

Target 

segment 

length 

(220MHz) 

Target 

segment 

length 

(300MHz) 

Target 

segment 

length 

(400MHz) 

1 

0.15A 

0.25A 

0.12A 

0.18A 

0.24A 

2 

0.15A 

0.25A 

0.12A 

0.18A 

0.24A 

3 

0.15A 

0.25A 

0.12A 

0.18A 

0.24A 

4 

0.15  A 

0.25A 

0.12A 

0.15A 

0.20A 

5 

0.15A 

0.25A 

O.lOA 

0.12A 

0.16A 

6 

0.15  A 

0.25A 

O.lOA 

O.lOA 

0.13A 

7 

0.15A 

0.25A 

0.12A 

0.16A 

0.21A 

8 

0.15A 

0.25A 

0.18A 

0.22A 

0.29A 

9 

0.15A 

0.25A 

0.1 8A 

0.22A 

0.29A 

10 

0.15A 

0.25A 

0.20A 

0.24A 

0.32A 

11 

0.15A 

0.25A 

0.20A 

0.24A 

0.32A 

12 

0.15A 

0.25A 

0.18A 

0.25A 

0.33A 

13 

O.lOA 

O.lOA 

0.16A 

0.22A 

0.29A 

O.lOA 

O.lOA 

O.lOA 

0.13A 

0.17A 

15 

O.lOA 

O.lOA 

O.lOA 

O.lOA 

0.13A 

16 

O.lOA 

O.lOA 

O.lOA 

0.15A 

0.2A 

17 

O.lOA 

O.lOA 

O.lOA 

0.15A 

0.2A 

Total 
No.  of 
segment 

1936 

1092 

4502 

4982 

4982 

Table  1;  The  segmentation  schemes  used  to  generate 
the  grid  model  of  the  aircraft  at  120MHz,  220MHz, 
300MHz  and  400MHz. 

The  segmentation  schemes  used  to  generate  the  numerical 
models  of  the  aircraft  at  220MHz,  300MHz,  and  400MHz, 
and  for  which  theoretical  results  are  presented  in  this  paper, 
are  indicated  in  Table  1.  The  total  number  of  segments  used 
to  numerically  model  the  aircraft  in  its  entirety  (as  indicated 
in  Table  1)  are  much  reduced  compared  to  that  which  would 
be  required  if  target  segments  length  of  0.1  A  were  used. 
Specifying  a  target  segment  length  of  0.1  A  would  result  in 
approximately  1 1400  segments  for  the  220MHz  grid  model, 
20000  segments  for  the  300MHz  grid  model  and  35300 
segments  for  the  400MHz  grid  model. 

4  PREDICTED  AND  MEASURED  RESULTS 

Appreciating  the  electrical  size  and  geometrical  complexity 
of  the  aircraft  model,  the  numerically  predicted  radiation 
patterns  are  compared  with  the  measured  radiation  patterns 
in  Fig.  6  to  Fig.  23.  A  close  examination  of  the  measured 
results  reveal  that  while  the  aircraft  is  symmetrical,  the 
measured  patterns  shows  some  asymmetry.  In  addition,  at 
angles  where  two  measurements  were  taken,  repeatability 
errors  are  noted;  the  discussion  in  Section  2  gave  the  reasons 
for  these  asymmetry  and  repeatability  errors.  Mostly  these 
repeatability  errors  are  noted  at  points  where  the  model  was 
tilted  backwards  and  the  positioner  interaction  as  well  as 
some  bending  of  the  perspex  mount  caused  measurement 
errors.  Although  comparison  of  patterns  in  terms  of 
normalized  values,  as  presented  in  [3,  4  and  5],  are  based 
on  pattern  integration  on  a  complete  volumetric  data  set, 
including  both  polarizations  for  both  measured  and  computed 
patterns,  the  comparison  between  measured  and  theoretical 
results  in  this  paper  is  shown  in  terms  of  absolute  gains 
which  are  associated  with  the  dominant  polarization.  The 
predicted  results  obtained  for  the  lower  frequencies 
numerical  models  exhibit  better  qualitative  agreement  with 
measured  results  because: 

•  electrically  smaller  size  of  the  overall  problem  tends  to 
result  in  simpler  pattern  shapes  with  fewer  main  lobes 
and  nulls. 

•  they  differ  in  their  segmentation  schemes  (as  shown  in 
Table  1)  and  generally  have  segmentation  closer  to  that 
which  is  required  by  the  NEC2  modelling  guide-lines 
[1]. 

•  they  are  generally  speaking  less  prone  to  numerical 
errors  due  to  the  smaller  problem  size. 
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Hie  400  MHz  cases  (Figs.  18  to  23)  shows  some  cases  where 
agreemmt  between  measurements  and  NEC2  results  are 
poor,  mainly  due  to  extreme  under-segmentation  (0.33X  for 
some  portions)  and  probably  represents  the  limit  to  which 
this  technique  can  be  stretched.  For  the  purpose  of  statistical 
link  analysis  these  results  were  still  useful,  but  they  may  be 
unsuitable  for  some  more  demanding  evaluations. 

5  CONCLUSION 

The  main  conclusions  and  recommendations  of  this  stu^  are 
listed  below: 

•  Under-segmentation  of  the  aircraft  geometry  in  areas 
removed  from  the  anteimas  made  the  evaluation  of 
antennas  at  UHF  frequencies  possible;  without  this 
technique  most  available  computer  resources  will  be 
inadequate  for  such  an  evaluation. 

•  Limited  measurements  (with  some  significant  errors  in 
some  cases)  proved  useful  in  qualitatively  assessing  the 
merit  of  the  computed  results. 

•  Measured  results  can  be  considerably  improved  by 
changing  the  mounting  method  on  scale  modeb  such 
that  only  azimuth  rotation  is  required.  This  can  be 
achieved  by  providing  side  mounts  as  well  as  front  and 
back  mounts  for  pitch  and  roll  plane  measurements 
respectively.  This  technique  was  used  during  subsequent 
smdies  with  different  aircraft  and  eliminated  the 
problem  of  repeatable  results  for  the  same  measurements 
points  (at  (f  and  360°,  for  instance). 

•  The  SIG  program  [2]  proved  to  be  very  useful  for 
automatic  grid  generation  of  the  aircraft  Particularly 
due  to  the  ability  to  vary  grid  size,  handle  symmetry 
and  generate  grid  models  for  different  frequencies. 
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Figure  1;  Hieoretical  azimuth-plane  radiation 
patterns  of  the  Top-Fin  antenna  at  120MHz  as 
obtained  by  NEC2  for  two  different  grid  models  of 
the  fighter  aircraft  Both  numerical  models  use  O.lOX 
segment  lengths  in  the  vicinity  of  the  Top-Fin 
antenna.  In  the  regions  which  are  sufficiently 
removed  from  the  Top-Fin  antenna,  the  one  model 
uses  0.15X  segment  lengths  while  the  other  model 
uses  0.25X,  segment  lengths. 
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Figure  2;  Hieoretical  pitch  elevation-plane  radiation 
patterns  of  the  Top-Fin  antenna  at  120MHz  as 
obtained  by  NEC2  for  two  different  grid  modeb  of 
the  fighter  aircraft  Both  numerical  modeb  use  O.lOX, 
segment  lengths  in  the  vicinity  of  the  Top-Fin 
antenna.  In  the  regions  which  are  sufficiently 
removed  from  the  Top-Fin  antenna,  the  one  model 
uses  0.15 A,  segment  lengths  while  the  other  model 
uses  0.25A,  segment  lengths. 


Figure  3:  The  user  defined  cross-sections  used  to 
generate  the  grid  model  of  the  fighter  aircraft  (The 
longitudinal  lines  shown  in  tliis  figure  are  just  for 
ease  of  vbualizing  the  structure  and  to  give  an 
indication  of  the  linear  interpolation  between  the 
cross  sections). 
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Figure  4;  The  grid  model  of  the  fighter  aircraft  at 
220MHz.  Since  symmetiy  is  used  for  simulations, 
only  half  of  the  grid  model  is  dbplayed.  The  two 
halves  are  connected  at  the  aircraft  axis  of  symmetry 
where  segments  are  seemed  to  be  terminated.  This 
numerical  model  comprises  of  4502  segments  with 
the  shortest  segment  lengths  being  O.lOA,  (in  the 
vicinity  of  the  antennas)  and  the  longest  segment 
lengths  being  0.20X  (in  a  region  removed  from  the 
antennas). 


Figure  5:  The  grid  model  of  the  fighter  aircraft  at 
300MHz.  Since  symmetry  is  used  for  simulations, 
only  half  of  the  grid  model  is  displayed.  The  two 
halves  are  coimected  at  the  aircraft  axis  of  symmetry 
where  segments  are  seemed  to  be  terminated.  This 
numerical  model  comprises  of  4982  segments  with 
the  shortest  segment  lengths  being  O.lOA,  (in  the 
vicinity  of  the  anteimas)  and  the  longest  segment 
lengths  being  0.25X,  (in  a  region  removed  from  the 
antennas). 
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Figure  6;  Measured  and  theoretical  azimuth 
radiation  patterns  of  the  Bottom-Fin  antenna  at 
220MHz.  The  numerical  model  comprises  of  4502 
segments  with  the  shortest  and  longest  segment 
lengths  being  O.lOX  and  0.20X  respectively. 


Figure  9;  Measured  and  theoretical  azimuth 
radiation  patterns  of  the  Top-Fin  antenna  at  220MHz. 
The  numerical  model  comprises  of  4502  segments 
with  the  shortest  and  longest  segment  lengths  being 
0.1  OA,  and  0.20 A,  respectively. 


Figure  7;  Measured  and  theoretical  side  roll 
elevation-plane  radiation  patterns  of  the  Bottom-Fin 
antenna  at  220MHz.  The  numerical  model  comprises 
of  4502  segments  with  the  shortest  and  longest 
segment  lengths  being  0.10^  and  0.20A.  respectively. 


Figure  10;  Measured  and  theoretical  side  roll 
elevation-plane  radiation  patterns  of  the  Top-Fin 
antenna  at  220MHz.  The  numerical  model  comprises 
of  4502  segments  with  the  shortest  and  longest 
segment  lengths  being  O.lOA,  and  0.20A,  respectively. 


Figure  8;  Measured  and  theoretical  pitch  roll 
elevation-plane  radiation  patterns  of  the  Bottom-Fin 
antenna  at  220MHz.  The  numerical  model  comprises 
of  4502  segments  with  the  shortest  and  longest 
segment  lengths  being  O.lOX  and  0.20A,  respectively. 


Figure  11;  Measured  and  theoretical  pitch  roll 
elevation-plane  radiation  patterns  of  the  Top-Fin 
antenna  at  220MHz.  The  numerical  model  comprises 
of  4502  segments  with  the  shortest  and  longest 
segment  lengths  being  0.10 A,  and  0.20 A  respectively. 
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Figure  12;  Measured  and  theoretical  azimuth 
radiation  patterns  of  the  Bottom-Fin  antenna  at 
300MHz.  The  numerical  model  comprises  of  4982 
segments  with  the  shortest  and  longest  segment 
lengths  being  O.lOX,  and  025^  respectively. 


Figure  15;  Measured  and  theoretical  azunuth 
radiation  patterns  of  the  Top-Fin  antenna  at  300MHz. 
The  numerical  model  comprises  of  4982  segments 
with  the  shortest  and  longest  segment  lengths  being 
O.lOA,  and  0.25X,  respectively. 
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Figure  13:  Measured  and  theoretical  side  roll 
elevation-plane  radiation  patterns  of  the  Bottom-Fin 
antenna  at  300MHz.  The  numerical  model  comprises 
of  4982  segments  with  the  shortest  and  longest 
segment  lengths  being  O.IOX.  and  0.25X.  respectively. 


Figure  16;  Measured  and  theoretical  side  roll 
elevation-plane  radiation  patterns  of  the  Top-Fin 
antenna  at  300MHz.  The  mnnerical  model  comprises 
of  4982  segments  with  the  shortest  and  longest 
segment  lengths  being  O.lOA,  and  0.2SA.  respectively. 
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Figure  14;  Measured  and  theoretical  pitch  roll 
elevation-plane  radiation  patterns  of  the  Bottom-Fin 
antenna  at  300MHz.  The  numerical  model  comprises 
of  4982  segments  with  the  shortest  and  longest 
segment  lengths  being  O.lOX,  and  0.25X  respectively. 


Figure  17;  Measured  and  theoretical  pitch  roll 
elevation-plane  radiation  patterns  of  the  Top-Fin 
antenna  at  300MHz.  The  numerical  model  comprises 
of  4982  segments  with  the  shortest  and  longest 
segment  lengths  being  O.lOA,  and  0.25A,  respectively. 
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Figure  18;  Measured  and  theoretical  azimuth 
radiation  patterns  of  the  Bottom-Fin  antenna  at 
400MHz.  Hie  numerical  model  comprises  of  4982 
segments  with  the  shortest  and  longest  segment 
lengths  being  0.13X  and  0.33X  respectively. 


Figure  21;  Measured  and  theoretical  azimuth 
radiation  patterns  of  the  Top-Fin  anteima  at  400MHz. 
The  numerical  model  comprises  of  4982  segments 
with  the  shortest  and  longest  segment  lengths  being 
0.13A,  and  0.33^  respectively. 


Figure  19;  Measured  and  theoretical  side  roll 
elevation-plane  radiation  patterns  of  the  Bottom-Fin 
anteima  at  400MHz.  Tbe  numerical  model  comprises 
of  4982  segments  with  the  shortest  and  longest 
segment  lengths  being  0.13A,  and  0.33X,  respectively. 


Figure  22;  Measured  and  theoretical  side  roll 
elevation-plane  radiation  patterns  of  the  Top-Fin 
anteima  at  400MHz.  The  numerical  model  comprises 
of  4982  segments  with  the  shortest  and  longest 
segment  lengths  being  0.13A,  and  0.33X,  respectively. 


Figure  20:  Measured  and  theoretical  pitch  roll 
elevation-plane  radiation  patterns  of  the  Bottom-Fin 
antenna  at  400MHz.  The  numerical  model  comprises 
of  4982  segments  with  the  shortest  and  longest 
segment  lengths  being  0.13^  and  0.33A,  respectively. 


Figure  23:  Measured  and  theoretical  pitch  roll 
elevation-plane  radiation  patterns  of  the  Top-Fin 
antenna  at  400MHz.  The  numerical  model  comprises 
of  4982  segments  with  the  shortest  and  longest 
segment  lengths  being  0.13A,  and  0.33X  respectively. 
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Abstract  -  A  new  surface  integral  equation  formulation 
based  on  radiation  currents  is  presented.  Numerical 
problems  have  been  observed  when  treating  low  contrast 
scatterers  with  surface  integral  equation  formulations 
based  on  the  total  equivalent  currents.  These  problems 
result  because  the  total  equivalent  currents  are  large 
compared  with  the  radiation  currents  for  low  contrast 
materials.  The  surface  integral  equation  formulation 
presented  here  avoids  this  problem  by  solving  directly  for 
the  radiation  currents. 

L  INTRODUCTION 

This  paper  presents  a  new  surface  integral 
equation  formulation  for  determining  the  electromagnetic 
field  scattered  by  objects  with  low  contrast.  Low  contrast 
materials  are  characterized  by  a  relative  permittivity  and 
relative  permeability  close  to  unity.  Low  contrast  materials 
include  compressed  gases  and  commonly  used  plastic 
foams. 

When  used  in  conjunction  with  a  numerical 
integral  equation  solution,  the  conventional  surface  integral 
equation  formulations  can  result  in  a  poor  estimate  of  the 
field  scattered  by  a  low  contrast  object.  The  poor 
performance  of  the  standard  formulations  is  due  to  their 
being  cast  in  terms  of  the  total  equivalent  current  on  the 
boundary  siuface  of  the  scattering  body.  Only  a  portion  of 
the  total  equivalent  current  gives  rise  to  the  scattered  field. 
When  this  portion  is  small  (as  in  the  case  of  low  contrast 
scatterers),  errors  in  the  radar  cross  section  (RCS) 
calculation  can  occur.  The  new  integral  equation 
formulation  is  cast  only  in  terms  of  the  portion  of  the  total 
current  that  gives  rise  to  the  scattered  field;  hence,  the  new 
formulation  avoids  this  problem. 

Here  the  new  integral  equation  formulation  is 
implemented  in  a  moment  method  (MM)  program  for  two- 
dimensional  objects.  Results  of  the  numerical  method  are 
verified  by  comparison  with  series  solution  results  for 
circular  cylindrical  objects. 


n.  THE  SURFACE  INTEGRAL  EQUATIONS 

Consider  a  closed  dielectric  scatterer  in  free  space 
with  surface  A  and  illuminated  by  source  currents  in  the 
free  space  region.  For  a  common  choice  of  the  combination 
constants  a  set  of  combined  field  integral  equations  [1,2] 
for  this  scatterer  may  be  written  succinctly  as 

-[el  (  J)  +  EL  (M)];j^  -  n  X  [h^  (J)  +  (M)]  - 

EL(J‘,M‘)/q„+nxH^(T,M‘)  (1) 

just  inside  A 

-[el  (  J)  +  el  (M)]^  +  n  X  [h-  (  J)  +  H-  (M)] 

=  0.  (2) 

just  outside  A 

The  following  definitions  are  made  with  respect  to  (1)  and 

(2): 

A  -  Dielectric/free  space  boundary 

J,M  -  Equivalent  electric  and  magnetic 

surface  currents  on  A 

J',M’  -  Source  currents  of  the  incident  wave 
n  -  Outward  surface  normal  on  A 
E*,H"^-  Electric  and  magnetic  fields  calculated  in  an 
infinite  homogeneous  free  space  region 
E',H”-  Electric  and  magnetic  fields  calculated  in  an 
infinite  homogeneous  region  with  constitutive 
parameters  of  the  dielectric 
q„  -  Characteristic  impedance  of  free  space 

The  tangential  portion  of  the  electric  field  on  A  is 
determined  using  the  expression 

E^=-nxnxE  (3) 

Equation  (1)  corresponds  to  the  exterior  equivalent 
situation  where  the  scatterer  is  replaced  with  free  space  and 
the  equivalent  surface  currents  radiate  the  scattered  field 
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outside  the  boimdary  of  the  scatterer  and  the  negative  of  the 
incident  fields  inside  the  scatterer.  Equation  (2) 
corresponds  to  the  interior  equivalent  situation  where  the 
free  space  region  is  replaced  with  dielectric  and  the 
negative  of  the  equivalent  surface  currents  radiate  zero 
fields  in  the  exterior  region  and  the  same  fields  as  in  the 
original  scattering  problem  in  the  interior  region. 

The  total  equivalent  surface  currents  can  be 
partitioned  into  scattering  (radiation)  currents  and  incident 
currents  as  follows: 


M  =  M'  +M° 

where 

J°  =nxH"(j‘,M‘)  (5) 

and 

M"  =-nxE^(j‘,M‘)  .  (6) 

The  incident  currents  J^andM”  are  the  total  equivalent 
surface  currents  for  a  free  space  or  "phantom"  dielectric 
scatterer.  Thus,  the  incident  currents  must  radiate  zero 
scattered  field  in  the  exterior  region  of  the  exterior 
equivalent  problem.  Clearly,  for  an  actual  dielectric 
scatterer,  the  scattering  currents  alone  produce  the  total 
scattered  field  in  the  exterior  equivalent  problem. 

Substituting  (4)  into  (1)  and  making  use  of  the 
linearity  of  the  electric  and  magnetic  field  operators  to 
manipulating  the  result  yields: 

-[el.  ( J- )  +  E;„  (M' )]— -  n  X  [h  ^  ( J- )  +  H  ^  (M- )]  = 

T|o 

[EL.(JMVI')  +  EL.(J”,M")]^  (7) 

+nx[H"(J-,M‘)  +  H"(J",M”)] 
just  inside  A 

Because  the  "incident  currents"  produce  the  negative  of  the 
incident  fields  in  the  interior  region  of  the  exterior 
equivalent  problem  (7)  can  be  reduced  to 

-[eL.(J-)+EL(M')]— -fix[H^(J')  +  H^(M')] 

=  0.  (8) 

just  inside  A 


Equation  (8)  is  a  statement  of  the  fact  that  the  scattering 
currents  must  radiate  zero  fields  in  the  interior  of  the 
exterior  scattering  problem.  Substituting  (4)  into  (2)  yields: 

-[EL.(J')  +  EL.(M')]-!-  +  nx[H-(J‘)  +  H-(M')]  = 

[e^  ( J" )  +  E^  (M" )]  -!-  -  fi  X  [h-  ( J" )  +  H-  (M“ )]  (9) 

just  outside  A 

Equations  (8)  and  (9)  are  a  set  of  coupled  integral 
equations  which  can  be  solved  for  the  scattering  ciurents 
directly.  In  this  new  integral  equation  formulation  the 
excitation  is  moved  from  the  exterior  region  to  the  interior 
region.  The  excitation  in  the  new  formulation  is  dependent 
on  the  material  parameters  of  the  scatterer.  For  a 
"phantom"  dielectric  scatterer,  the  right  hand  side  of  (9) 
(the  excitation)  is  zero  so  that  the  scattering  currents  are 
identically  zero. 

The  radiation  current  integral  equation 
formulation  presented  here  is  developed  from  the 
conventional  combined  field  integral  equation  (CFIE) 
formulation.  A  similar  radiation  current  integral  equation 
formulation  can  also  be  obtained  for  the  PMCHW  [2,3] 
formulation  by  following  the  same  steps.  The  radiation 
current  integral  equation  formulation  is  also  easily 
extended  to  the  multiple  dielectric  scatterer  case. 

m.  NUMEMCAL  RESULTS 

To  demonstrate  the  radiation  current  integral 
equation  formulation  it  was  implemented  in  a  two- 
dimensional  moment  method  program  [4].  For  this 
program  (8)  and  (9)  were  specialized  to  the  transverse 
magnetic  (TM)  and  transverse  electric  (TE)  polarization 
cases.  The  2-D  MM  program  used  a  piecewise  linear 
approximation  to  the  contour  of  the  scatterer  and  pulse 
current  expansion  with  point  matching.  Segments  of  the 
piecewise  linear  representation  of  the  scatterer  are  referred 
to  as  zones. 

Because  the  program  originally  used  the 
conventional  CFIE  there  was  no  need  to  change  the 
calculation  of  the  MM  impedance  matrix  when 
implementing  the  new  formulation  (note  that  the  operator 
forms  on  the  left  sides  of  (1)  and  (2)  are  identical  to  those 
on  the  left  sides  of  (8)  and  (9)).  Only  the  calculation  of  the 
excitation  required  modification.  Accurate  calculation  of 
the  excitation  is  critical  to  the  success  of  the  method.  On 
linear  zones  the  incident  currents  J°  andM”  have  constant 
magnitude  and  linear  phase  for  an  incident  plane  wave.  If 
the  excitation  in  (9)  is  calculated  using  the  existing 
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impedance  elements  (which  assume  constant  magnitude 
and  constant  phase)  it  was  observed  that  the  radiation 
current  based  formulation  yields  no  better  results  than  the 
conventional  formulation  when  treating  low  contrast 
scatterers.  This  observation  is  true  for  both  the  TE  and  TM 
cases.  To  achieve  accurate  results,  the  phase  variation  of 
the  incident  currents  on  each  zone  must  be  incorporated  in 
the  calculation  of  the  excitation.  The  best  approach  is  to 
derive  new  impedance  elements  which  incorporate  the 
linear  phase  variation  of  the  incident  currents.  An 
alternative  numerical  approximation  is  to  subdivide  each 
linear  zone  into  an  odd  number  of  subzones  to  more 
accurately  model  the  phase  variation  of  the  incident 
currents.  In  this  approach  the  field  at  the  match  point  of  the 
center  subzone  of  each  zone  is  calculated  using  the  existing 
impedance  elements  (not  the  ones  in  the  MM  matrix  but 
ones  for  each  subzone  calculated  with  the  same  routines). 
The  use  of  an  odd  number  of  subzones  is  required  so  that 
match  point  of  the  center  subzone  will  be  the  same  as  the 
match  point  of  the  entire  zone.  The  former  approach  is 
more  accurate  and  requires  less  computation.  The  latter 
approach,  however,  has  the  advantage  of  not  being  limited 
to  plane  wave  excitation.  For  expedience,  the  latter 
approach  was  taken  here  with  each  zone  being  divided  into 
eleven  subzones  for  the  calculation  of  the  excitation. 

Figure  1  compares  the  bistatic  RCS  of  a  low 
contrast  dielectric  cylinder  calculated  using  the 
conventional  formulation  with  those  calculated  using  the 
radiation  current  formulation.  Also  plotted  in  Fig.  1  are 
series  solution  results.  The  radius  of  the  cylinder  is  one  free 
space  wavelength,  A,„.  The  dielectric  constant  of  the 
cylinder  is  1.01.  The  MM  program  is  written  in  FORTRAN 
and  uses  single  precision  floating  point  numbers.  For  both 
numerical  solutions  the  cylinder  was  divided  into  100 
linear  zones.  The  plane  wave  was  incident  from<t>  =  180" . 
In  Fig.  1  the  RCS,  a,  is  normalized  by  .  Figure  la  is  for 
a  TM  polarized  incident  wave  and  Fig.  lb  is  for  a  TE 
polarized  incident  wave.  Figure  1  demonstrates  the  error  in 
the  calculated  RCS  which  can  occur  when  treating  low 
contrast  scatterers  using  the  conventional  formulation. 

Figure  2  compares  the  magnitude  of  the 
normalized  electric  and  magnetic  scattering  current  for  the 
case  of  Fig.  lb  (TE  polarization).  Series  solution  results  are 
also  given.  Figure  2a  plots  the  electric  current  and  Fig.  2b 
plots  the  magnetic  cinrent.  The  electric  current  is 
normalized  by  the  incident  magnetic  field  and  the  magnetic 
current  is  normalized  by  the  incident  electric  field.  Because 
the  conventional  formulation  calculates  the  total  equivalent 
currents,  the  scattered  currents  in  Fig.  2  for  the 
conventional  formulation  were  determined  by  subtracting 
the  incident  currents  from  the  calculated  total  currents.  The 
radiation  current  formulation  determines  the  scattered 


currents  directly.  Figure  2  illustrates  that  use  of  the 
conventional  formulation  results  in  errors  in  the  scattering 
currents  for  low  contrast  scatterers.  As  expected  these 
errors  occur  at  places  on  the  scatterer  where  the  scattering 
current  is  small. 

IV.  CONCLUSIONS 

A  radiation  current  based  surface  integral  equation 
formulation  for  low  contrast  scatterers  was  presented. 
Numerical  results  demonstrate  the  utility  and  validity  of  the 
new  formulation. 

Although  its  main  application  appears  to  be  in 
treating  low  contrast  scatterers,  the  radiation  current 
formulation  is  also  useful  for  accurately  determining  the 
equivalent  currents  on  general  material  bodies  in  regions  of 
the  siuface  where  the  radiation  currents  are  small.  Note 
that  unlike  the  most  commonly  used  low-contrast 
formulation,  the  Muller  formulation  [3],  the  radiation 
current  formulation  works  for  high  contrast  scatterers  as 
well  as  low  contrast  scatterers. 

Implementation  of  the  formulation  in  existing 
moment  method  programs  is  relatively  straightforward 
since  the  moment  matrix  is  unchanged.  Only  the  excitation 
vector  must  be  changed. 

REFERENCES 

[1]  J.  R.  Mautz  and  R.  F.  Harrington,  "Boimdary 
formulations  for  aperture  coupling  problems," 
Arch.  Elek.  Ubertragung,  vol.  34,  pp.  377-384, 
1980. 

[2]  A.  A.  Kishk  and  L.  Shafai,  "On  the  accuracy 
limits  of  different  integral  equation  formulations 
for  numerical  solution  of  dielectric  bodies  of 
revolution,"  Can.  J.  Phys.,  vol.  63,  pp.  1532-1539, 
1985. 

[3]  J.  R.  Mautz  and  R.  F.  Harrington, 
"Electromagnetic  scattering  from  a  homogeneous 
body  of  revolution,"  Arch.  Elek.  Ubertragung,  vol. 
32,  pp.  71-80,  1979. 

[4]  P.  M.  Goggans  and  T.  H.  Shumpert,  "CFIE  MOM 
solution  for  TE  and  TM  incidence  on  a  2-D 
conducting  body  with  a  dielectric  filled  cavity," 
IEEE  Trans.  Antennas  Propagat,  vol.  AP-38,  pp. 
1645-1649,  October  1990. 


17 


(a)  TM  polarization 


(b)  TE  polarization 


Figure  1.  Bistatic'  RCS  of  a  dielectric  cylinder  with  a 
radius  of  one  ffee-space  wavelength  and  a 
dielectric  constant  of  1.01.  The  plane  wave  is 
incident  from  <|)  =  180° . 


(a)  Electric  current 


(b)  Magnetic  current 


Figure  2.  Normalized  scattering  currents  on  of  a 
dielectric  cylinder  with  a  radius  of  one  free- 
space  wavelength  and  a  dielectric  constant  of 
1.01.  The  TE  polarized  plane  wave  is  incident 
from  (j)  =  180° . 
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The  Measured  Equation  of  Invariance  Method  Applied  to  Randomly 

Rough  Surfaces 

John  B.  Schneider  and  Shira  Lynn  Broschat* 


Abstract 

The  Measured  Equation  of  Invariance  (MEI)  method  has 
received  considerable  attention  recently.  Unlike  more 
traditional  numerical  techniques,  MEI  solutions  cire  ob- 
tcdned  by  inversion  of  a  relatively  small  and  sparse  matrix. 
Therefore,  the  MEI  method  can  potentially  provide  a  so¬ 
lution  much  more  quickly  than  other  techniques.  To  date, 
the  MEI  method  has  been  applied  primarily  to  discrete 
objects.  In  this  paper,  bistatic  radar  cross  sections  for 
one-dimensional,  perfectly  conducting,  randomly  rough 
smfaces  are  obtained  using  the  MEI  method.  The  imple¬ 
mentation  suitable  for  this  problem  requires  some  mod¬ 
ification  and  enhancement  of  the  original  algorithm  to 
achieve  the  desired  accuracy.  These  algorithmic  changes 
can  be  applied  to  the  discrete  scattering  problem  as  well. 
Monte  Carlo  results  for  the  bistatic  scattering  cross  sec¬ 
tion  for  surfaces  with  Gaussian  statistics  and  satisfying  a 
Gaussian  roughness  spectrum  are  compared  to  those  from 
another  technique  and  excellent  agreement  is  obtained. 

I.  Introduction 

The  Measured  Equation  of  Invariance  (MEI)  method 
was  recently  introduced  [1,  2,  3,  4]  cis  a  way  to  deter¬ 
mine  the  electromagnetic  fields  scattered  from  discrete 
objects.  Initially,  results  were  reported  mostly  for  con¬ 
ducting  two-dimensional  objects,  but  the  method  can  be 
applied  much  more  broadly.  Hybrid  techniques  have  been 
used  in  conjunction  with  the  MEI  method  to  obtain  scat¬ 
tering  from  penetrable  objects  [5,  6,  7,  8];  however,  many 
of  these  scatterers  can  be  analyzed  using  a  simpler  scheme 
if  the  metrons  are  selected  carefully  [2].  Additionally,  the 

‘The  authors  are  with  the  School  of  Electrical  Engineering  and 
Computer  Science,  Washington  State  University,  Pullman,  WA 
99164-2752 


MEI  method  can  be  used  to  provide  solutions  to  Laplace’s 
equation  [3,  9]. 

The  principal  attraction  of  the  MEI  method  is  that  it 
yields  a  solution  via  a  sparse  and  relatively  small  ma¬ 
trix,  and  hence  requires  much  less  time  and  memory  than 
more  conventioncil  techniques.  The  computationcJ  mesh 
is  terminated  very  close  to  the  scattering  object.  Interior 
mesh  points  are  related  to  each  other  using  standard  finite 
difference  techniques  while  boundary  nodes  are  hzindled 
using  the  Measured  Equations  of  Inveiriance.  The  time 
required  to  fill  the  matrix  is  roughly  O(A^),  where  N  is 
the  number  of  unknowns  on  the  boundary.  Because  the 
matrix  is  sparse,  the  time  required  to  invert  the  matrix  is 
small  {0{N)y,  thus,  the  total  computation  time  is  domi¬ 
nated  by  the  time  required  to  fill  the  matrix.  This  is  in 
contrast  to  the  traditional  Method  of  Moments  approach 
which  requires  the  inversion  of  a  full  matrix  (an  0{N^) 
operation). 

In  this  paper  the  MEI  method  is  used  to  solve  for  the 
bistatic  radar  cross  sections  for  randomly  rough  surfaces. 
Specifically,  ensemble  averaging  over  a  number  of  surface 
realizations  is  used  to  approximate  the  results  for  a  sin¬ 
gle,  infinitely  long  surface.  Surfaces  with  Gaussian  statis¬ 
tics  and  satisfying  a  Gaussian  roughness  spectrum  are 
assumed.  However,  it  should  be  emphasized  that  this  nu¬ 
merical  technique  is  not  restricted  to  this  class  of  surfaces. 

In  Section  II  a  brief  review  of  the  MEI  method  is  pro¬ 
vided.  The  pertinent  aspects  of  the  rough  surface  scatter¬ 
ing  problem  are  outlined  in  Section  III.  In  Section  IV  our 
implementation  of  a  MEI-based  solution  to  the  problem 
is  described.  Finally,  results  are  presented  in  Section  V. 

II.  Review  of  the  MEI  Method 

Ultimately,  our  goal  is  to  solve  the  wave  equation  sub¬ 
ject  to  the  appropriate  boundary  conditions.  For  the 
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sake  of  concreteness,  we  assume  TM  illumination  of  a 
perfectly-conducting,  two-dimensional  scatterer  (e.g.,  a 
surface  for  which  the  height  varies  as  a  function  of  one 
spatial  coordinate)  so  that  the  governing  differential  equa¬ 
tion  is  reduced  to  a  scalar  one: 

(v2  -1-  e)E,  =  0  (1) 

where  k  is  the  wavenumber  of  the  incident  field.  If  space  is 
discretized  into  a  uniform  mesh  and  a  local  configuration 
of  nodes  is  as  shown  in  Fig.  1,  (1)  can  be  approximated 
using  central  differences  to  obtain 

-  (4  -  +  E,i  +  E,2  +  E,3  +  E,4  =  0  (2) 

where  A  is  the  separation  between  mesh  points.  This  can 
be  written  more  generally  as 

N 

^  ^  —  0,  (3) 

i=0 

where  N  is  the  number  of  nodes  “connected”  to  the  ze¬ 
roth  (central)  node  and  the  aj’s  represent  appropriate 
weights.  For  a  Cartesian  structure,  such  as  in  Fig.  1,  all 
weights  for  the  surrounding  nodes  are  equal.  Had  a  po¬ 
lar  mesh  been  used,  the  weights  would  have  to  account 
for  the  global  location  within  the  mesh,  i.e.,  the  weights 
would  be  a  ftmction  of  the  radial  distance  from  the  ori¬ 
gin.  If  a  higher-order  differencing  scheme  had  been  used, 
N  would  have  to  be  increased.  These  difference  equations 
are  invariant  to  the  location  and  geometry  of  the  scatterer 
and  to  the  field  of  excitation.  The  appeal  of  the  finite  dif¬ 
ference  formulation  is  the  sparsity  and  ease  of  calculating 
the  non-zero  elements  in  the  resulting  matrix.  However, 
the  finite  difference  approach  is  problematic  in  that  it 
provides  no  simple  way  to  terminate  the  computational 
mesh  for  unbounded  problems. 

To  address  this  shortcoming,  Mei  developed  the  Mea¬ 
sured  Equation  of  Invariance  (MEI)  method  that  com¬ 
bines  features  of  both  differential  and  integral  based 
methods  [4].  Basically,  the  MEI  method  provides  a  means 
to  select  appropriate  a,  ’s  in  (3)  so  that  the  mesh  need  not 
be  orthogonal.  This  allows  the  fields  at  the  nodes  on  the 
edge  of  the  computational  domain  to  be  related  simply  to 
points  in  the  interior.  Furthermore,  Mei  maintains  that 
the  resulting  set  of  equations  (and  hence  the  a,’s)  is  lo¬ 
cation  dependent,  geometry  specific,  and  invariant  to  the 
field  of  excitation.  The  MEI  method  does  not  seek  to  find 


a  strict  discretization  of  the  wave  equation.  Instead,  the 
discretization  that  is  obtained  is  an  approximation  to  an 
unknown  operator  that  satisfies  both  the  wave  equation 
and  the  radiation  condition.  This  is  discussed  in  more 
detail  in  [2]. 

To  illustrate  the  technique,  consider  the  configuration 
of  nodes  shown  in  Fig.  2.  These  nodes  are  typical  of  the 
mesh  at  the  edge  of  the  computational  domain.  We  seek 
a,  ’s  such  that 

3 

^a.E„  =  0.  (4) 

*=0 

Since  (4)  is  a  homogeneous  equation,  one  of  the  weights 
(e.g.,  ao)  may  be  chosen  arbitrarily.  The  remaining  three 
weights  are  determined  via  three  equations.  These  equa¬ 
tions  are  obtained  by  assuming  independent  source  distri¬ 
butions,  known  as  metrons,  over  the  surface  of  the  scat¬ 
terer.  Each  metron  gives  rise  to  a  field,  known  as  the 
measured  field,  which  is  easily  calculated.  The  weights 
are  obtained  by  satisfying  (4)  for  each  of  the  measured 
fields.  Specifically,  three  metrons  are  assumed  J],  J^, 
and  Jg.  From  these  the  nth  measured  field  E"(t)  is  ob¬ 
tained  via 

E”(r)  =  f  J^{T')Gir\r')ds'  (5) 

J  s 

where  r  is  the  observation  point,  r'  is  a  point  on  the  sur¬ 
face,  and  G(r|r')  is  the  appropriate  Green’s  function  (typ¬ 
ically  the  free  space  Green’s  function).  Assuming  ao  =  1, 
the  remaining  weights  are  obtained  using 

■  ai  1  [  El,  El,  El,  1  r  Elo  ' 

a2  =-  El,  El,  El,  .  El,  (6) 

.  «3  J  El,  El,  El,  J  El,  _ 

which  Ccm  be  written  as 

a  =  -M“‘  •  Mo  (7) 

where  the  matrix  M  contains  the  measured  fields  at  the 
neighboring  nodes  while  the  vector  Mo  contains  the  fields 
sampled  at  the  zeroth  node. 

Although  the  metrons  are  chosen  to  be  independent, 
the  sampled  values  of  the  measured  fields  over  a  small 
number  of  mesh  points  may  not  be  independent.  Hence, 
the  matrix  M  may  be  singular  and  a  solution  to  (7)  will 
not  exist.  There  are  several  ways  to  circumvent  this  prob¬ 
lem.  Perhaps  the  easiest  is  to  use  more  metrons,  and  sub¬ 
sequently  more  measured  fields,  than  unknowns.  In  this 
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Figure  1;  Nodes  in  the  computational  mesh. 
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Figure  2:  Typical  configuration  of  nodes  at  the  edge  of  the  computational  mesh 


case,  the  weights  are  obtained  via  least  squares  so  that  a 
is  the  vector  that  minimizes  ||M  •  a  +  Mo||.  This  least 
squares  operation  is  performed  on  a  small  matrix  (e.g., 
in  the  work  presented  here  M  is  5  x  3)  and  thus  requires 
neglible  time  relative  to  the  calculation  of  the  measured 
fields. 

In  principle  the  MEI  method  can  be  used  for  cJl  nodes 
within  the  computational  domain.  However,  this  would 
require  careful  treatment  of  the  singularity  in  the  inte¬ 
grand  of  (5)  for  nodes  adjacent  to  the  surface  (i.e.,  at 
least  one  neighboring  node  would  be  on  the  surface  and 
the  measured  field  for  this  node  would  require  careful  cal¬ 
culation  of  (5)).  To  obviate  the  consideration  of  this  sin¬ 
gularity,  a  layer  of  nodes  that  use  standard  finite  differ¬ 
ence  coefficients  is  placed  adjacent  to  the  surface.  If  the 
mesh  is  non-orthogonal,  but  relatively  undistorted,  Pous 
has  found  that  weights  can  be  obtcdned  using  a  local  po¬ 
lar  approximation  [2].  This  layer  of  nodes  insures  that 
the  measured  fields  only  need  to  be  calculated  at  points 
above  the  surface. 

Prom  this  point,  the  MEI  method  is  similar  to  the  fi¬ 
nite  difference  method.  The  weights  found  above  are  used 
to  construct  a  “connectivity  matrix”  A  that  specifies  the 
relationship  between  every  node  in  the  mesh.  Each  row 
number  corresponds  to  the  global  node  number  for  a  given 
unknown  (i.e.,  the  field  at  that  node).  The  non-zero  ele¬ 
ments  in  that  row  correspond  to  the  a.’s  associated  with 


the  neighboring  nodes  cind  the  node  itself.  The  connec¬ 
tivity  matrix  is  sparse,  and  it  may  have  other  structure 
(such  as  bandedness)  that  can  be  exploited.  For  nodes 
that  are  adjacent  to  the  surface  of  the  scatterer,  one  of 
the  neighboring  values  corresponds  to  the  known  surface 
field  which  serves  as  the  forcing  function  to  the  system  of 
equations.  Symbolically,  this  is  written  as 

A  •  E  =  F  (8) 

where  E  contains  the  unknown  field  values  at  each  node 
in  the  computational  mesh  and  F  is  the  forcing  function 
determined  by  a  combination  of  the  known  surface  field 
and  the  weights  determined  using  the  loccil  polai’  approx¬ 
imation.  The  solution  is  then  obtained  via 

E  =  A-*F.  (9) 

Once  E  has  been  determined,  the  actual  surface  currents 
are  deduced  and  far-field  quantities  are  obtained  as  de¬ 
scribed  in  Sec.  IV. 

III.  Rough  Surface  Scattering 
Problem 

Currently  rough  surface  scattering  is  of  interest  to  re¬ 
searchers  in  a  variety  of  disciplines.  It  has  applications 
in  such  diverse  areas  as  ultrasonics,  radar  imaging,  sonar 
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detection,  solid-state  physics,  optics,  astronomy,  and  mi¬ 
crowave  remote  sensing. 

In  the  past  three  decades,  much  work  has  been  done 
in  the  development  of  approximate  analytic  models  to 
predict  wave  scattering  from  rough  surfaces.  These  in¬ 
clude  the  small  slope  approximation  [10,  11,  12],  the 
phase  perturbation  technique  [13,  14],  the  operator  ex¬ 
pansion  method  [15,  16],  the  Dashen-Wurmser  approx¬ 
imation  [17],  the  unified  perturbation  method  [18,  19], 
and  the  quasi-slope  approximation  [20].  In  addition,  some 
work  has  been  done  in  the  development  of  Monte  Carlo 
numerical  techniques  that  are  exact  in  the  sense  that  no 
physical  approximations  are  made  in  the  underlying  equa¬ 
tions.  These  include  techniques  based  on  integral  equa¬ 
tions  [21,  22],  the  finite  element  method  [23],  and  finite 
difference  methods  [24,  25,  26].  In  this  work  we  consider 
the  MEI  method.  We  restrict  our  consideration  to  per¬ 
fectly  conducting,  one-dimensional,  randomly  rough  sur¬ 
faces  with  Gaussian  statistics  eind  satisfying  a  Gaussian 
surface  roughness  spectrum  [27].  Surfaces  with  Gaussian 
statistics  have  been  studied  extensively  [11,  21,  28,  29]. 

The  problem  geometry  is  shown  in  Fig.  3.  TM  illumi¬ 
nation  is  assumed  so  the  total  electric  field  at  the  surface 
vanishes  (Dirichlet  boundary  condition).  The  surface  pro¬ 
file  is  given  by  f{x).  The  mean  height  of  the  surface  is 
zero,  i.e.,  {f)^  =  0,  where  {•)^  indicates  averaging  over  the 
entire  surface.  The  standard  deviation,  or  RMS  surface 
height,  is  given  by  /i  =  y/{P)~-  3°*^  surface  heights  and 
slopes  are  Gaussian  distributed. 

The  normalized  correlation  function  is  defined  by 

C(c)  =  + 

which  for  Gaussian  surfaces  is  given  by 

C(C)  =  exp  . 

The  correlation  length  I  is  the  length  at  which  the  cor¬ 
relation  fimction  decreases  by  a  factor  of  1/e.  The  sur¬ 
face  roughness  spectrum  is  obtained  by  taking  the  Fourier 
transform  of  h^C(Q.  For  Gaussian  surfaces,  the  statis¬ 
tics  are  completely  specified  by  just  two  parameters,  the 
correlation  length  and  the  RMS  surface  height,  and  gen¬ 
eration  of  siuface  realizations  is  relatively  simple  [21,  27]. 

In  numerical  simulations  of  rough  surface  scattering, 
finite-length  surfaces  must  be  used  to  model  scattering 
from  infirute  surfaces.  When  a  single  plane  wave  strikes 


(10) 

(11) 


a  finite-length  surface,  edge  diffraction  occurs.  One  way 
of  minimizing  diffraction  effects  is  to  construct  an  inci¬ 
dent  wave  that  tapers  to  very  small  values  at  the  surface 
edges.  Diffraction  still  occurs,  but  it  makes  negligible  con¬ 
tributions  to  the  scattered  field.  Tapered  incident  waves 
have  been  introduced  by  Thorsos  [21]  and  Chan  and  Fung 
[22].  The  tapered  incident  field  used  by  Thorsos  is  an  ap¬ 
proximation  to  an  incident  field,  consisting  of  an  angular 
spectrum  of  plane  waves,  that  satisfies  the  wave  equation 
exactly.  This  approximate  field  is  given  by 

EI{t)  =  exp  {-jki  •  r  [1  -t-  u;(r)]  -  (a;  -ytan9i)/g^} 

(12) 

where 

u;(r)  =  [2(x  -  ytaxiOif/g^  -  l]  /{kgoosOi)^,  (13) 

r  =  {x,y)  is  a  point  above  the  surface,  0,-  is  the  in¬ 
cident  angle  measured  from  the  vertical,  and  k,-  = 
ko(sm  Oi,  cos 6i)  is  the  free-space  incident  wave  vector  in 
the  xy  plane.  Equation  (12)  satisfies  the  wave  equation 
to  order  l/(fcogcos0,)2  for  kog  cos  6i  »  1.  The  parame¬ 
ter  g  controls  the  tapering,  aind  care  must  be  taken  in  its 
choice.  Angular  resolution,  edge  effects,  and  accuracy  in 
satisfying  the  wave  equation  all  depend  on  g  [21].  In  ad¬ 
dition,  the  tapering  must  be  accomplished  in  such  a  way 
that  differences  between  the  finite  surface,  tapered  plane 
wave  results  and  infinite  surface,  single  plane  wave  results 
are  negligible.  For  the  numerical  examples  presented  in 
this  paper,  g  =  L/4  is  used,  where  L  is  the  horizontal 
extent  of  each  surface. 

For  each  numerical  study,  50  finite-length  surfaces  are 
generated  using  the  method  proposed  by  Thorsos  [21]. 
Monte  Carlo  results  are  then  obtained  by  taking  the  en¬ 
semble  average  of  the  cross  sections  of  the  50  surfaces. 
The  general  procedure  is  to  randomly  generate  a  surface 
spectrum  that  has  Gaussian  statistics  and  then  inverse 
transform  the  spectrum  to  obtain  a  surface  profile.  Each 
surface  consists  of  N  discrete  points  horizontally  sepa¬ 
rated  by  Aa;  over  a  surface  length  L'.  The  x  component 
of  each  point  along  the  surface  is  specified  by  the  location 
x„  =  nAx  for  1  <  n  <  iV.  The  surfaces  are  generated 
using 

j  N/2-\ 

/(®n)  =  77  XI  ^{I^e)^^Y>\-3Kex„\  (14) 

^  t--NI2 
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Figure  3:  Problem  geometry. 


where,  for  0  <  ^  <  N/2, 

F{Ke)  =  V2^L>W{Ke)^{N{0, 1)  -  jA(0, 1)]  (15) 

and  for  f  =  0  or  N/2 

F{Kt)  =  v/27rLW(/^£)A(0, 1)  (16) 

In  (14)-(16) 

h^l 

WiKe)  =  ^  exp[-/C|zV4]  (17) 

is  the  Gaussian  surface  roughness  spectrum,  Ke  = 
27^1/ L',  h  is  the  RMS  surface  height,  I  is  the  correlation 
length,  and  A(0, 1)  is  a  number  sampled  from  a  Gaus¬ 
sian  distribution  with  zero  mean  and  unity  variance.  For 
e  <  0,  FiK-e)  =  F{Ke)*. 

The  nature  of  the  discrete  spectrum  causes  correlation 
of  the  ends  of  each  surface.  To  circumvent  this,  an  ex¬ 
tended  surface,  several  times  longer  than  the  N  required 
points,  is  generated.  Each  surface  used  in  the  numeri¬ 
cal  simulation  is  then  cut  from  the  longer  surface  and, 
hence,  the  correlation  of  the  ends  is  negligible.  For  the 
numerical  studies  presented  here,  surfaces  with  a  length  of 
L'  =  256A  were  generated  but  segments  of  length  L  =  80A 
were  used  in  the  calculations  (A  is  the  wavelength  of  the 
illumination). 

IV.  MEI  Method  for  Rough 
Surfaces 

In  order  to  use  the  MEI  method,  a  computational  mesh 
enclosing  the  scattering  object  must  be  specified.  Figure  4 


shows  a  segment  of  a  typical  surface  with  a  two-layer  com¬ 
putational  mesh.  To  generate  the  mesh,  two  mesh  points 
are  specified  for  each  surface  point.  Both  points  are  along 
the  surface  normal  and  are  separated  by  a  distance  At). 
To  obtain  the  surface  normal,  the  smrface  slope  is  needed. 
The  X  and  y  components  of  the  surface  unit  normal  vector 
are 

n.  = 

1 

y/l  +  [f'ix)f 

where  f'{x)  =  df/dx.  The  surface  derivative  could  be 
obtained  approximately  using  finite  differences;  however, 
the  surface  spectrum  is  needed  to  generate  the  surface 
realizations  and  is  avciilable  to  obtain  surface  slopes  (see 
(14)).  Thus,  we  find  f'{x)  via 

nx)  =  T-\jKF{K))  (20) 

where  denotes  the  inverse  Fourier  trcinsform. 

Note  that  the  surface  generation  scheme  creates  points 
that  are  horizontally  offset  from  neighboring  points  by 
Ax.  However,  since  neighboring  points  do  not  necessarily 
have  the  Scime  vertical  components,  the  distance  between 
points  is  nonuniform. 

Figure  5  shows  cin  expanded  view  of  one  comer  of  the 
computationcil  mesh.  The  global  node  numbers  are  shown 
together  with  some  of  the  node  interconnections.  Odd- 
numbered  nodes  correspond  to  nodes  in  the  interior  while 
even-numbered  nodes  are  on  the  outer  boundary.  The 
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first  and  lEist  two  nodes  are  exceptional  and  must  be 
given  special  consideration.  In  our  implementation,  all 
odd  nodes  (except  1  and  iV  —  1)  are  connected  to  neigh¬ 
boring  nodes  in  a  manner  similar  to  node  13  shown  in 
Fig.  5.  Since  neighbors  are  available  above,  below,  to  the 
left,  and  to  the  right,  the  local  polar  scheme  proposed  by 
Pous  can  be  employed  [2].  As  mentioned  before,  this  sim¬ 
plifies  calculation  of  the  connectivity  weights  and  makes 
consideration  of  the  singularity  in  the  Green’s  function  at 
the  surface  unnecessary.  All  even  nodes  between  2  and 
N  are  handled  using  the  MEI  method.  These  nodes  are 
connected  to  one  neighbor  in  the  interior  and  to  the  neigh¬ 
bors  to  the  left  and  right  as  shown  for  node  8.  Nodes  2 
and  N  are  both  connected  to  only  two  neighbors  and  are 
handled  using  the  MEI  method.  Nodes  1  and  N  —  1  are 
treated  using  the  local  polar  approximation.  For  these 
nodes,  an  additional  node  is  assumed  to  exist  to  the  left 
of  node  1  or  to  the  right  of  node  N  —  1  so  that  the  proper 
weights  can  be  determined.  However,  the  field  is  assumed 
to  be  zero  at  the  fictitious  neighbor  and  it  does  not  enter 
into  the  calculation.  Note  that  the  tapered  incident  field 
is  small  near  nodes  1  and  AT  - 1  so  setting  the  field  to  zero 
should  introduce  negligible  error.  With  this  connectivity 
scheme,  the  connectivity  matrix  A  is  not  only  sparse,  but 
is  also  tightly  banded — there  are  a  total  of  five  diagonals 
(including  the  main  diagonal)  that  contain  non-zero  ele¬ 
ments.  A  number  of  public  domain  routines  exist,  such 
as  those  in  LINPACK  [30],  that  can  be  used  to  efficiently 
invert  such  a  matrix. 

In  many  rough  surface  scattering  problems  bistatic 
results  span  a  range  of  over  60  dB.  For  instance,  for 
slightly  rough  surfaces  the  coherent  specular  reflection 
contains  significcint  energy  relative  to  the  incoherent  en¬ 
ergy  found  at  grazing  angles.  This  dynamic  range  re¬ 
quires  an  extremely  accurate  numerical  solution.  In  pre¬ 
viously  published  results  the  accuracy  of  MEI-based  so¬ 
lutions  was,  at  most,  weakly  dependent  on  the  selection 
of  metrons.  However,  we  found  that  using  an  initial  se¬ 
lection  of  metrons  that  are  physically  unrealizable  leads 
to  poor  results.  Because  of  the  accuracy  required  Euid  the 
difficulty  in  making  a  “good”  initial  selection  of  metrons, 
an  iterative  scheme  was  used.  The  initial  set  of  metrons 

WclS 

Jj(r,)  =  expi-xf/g^) 

-  cos(fcoC.) 


J®(r,)  =  sin(feo^t) 

J^(ri)  =  exp(-a;?/fli^)cos(fcoyt) 

J®(r,)  =  exp(-x?/fl|2)sin(fcoyi) 

where  g  is  the  taper  factor  commensurate  with  the  inci¬ 
dent  field,  r;  =  (x,',y,  )  corresponds  to  the  ith  point  along 
the  surface,  is  the  path  length  along  the  surface  from 
the  first  point  to  the  ith  point,  and  fco  is  the  free-space 
wavenumber.  Using  the  following  Green’s  function 

G{r\r')  =  ^H^^\ko\v-T'\),  (21) 

where  rjo  is  the  characteristic  impedance  of  free  space,  the 
measured  fields  cire  obtained  (see  (5))  using 

N 

E^ir)  «  ^r(r.)G(r|r.)Ae.-  (22) 

i=l 

where  A^,-  is  the  distance  from  the  ith  to  the  (i  +  l)th 
point  (A^at  is  set  equal  to  A^Af-i)- 

The  iterative  scheme  proceeds  as  follows.  The  metrons 
shown  above  are  employed  in  the  MEI  method  to  obtain 
the  fields  above  the  surface.  From  these  fields  a  surface 
current  is  obtained.  The  calculated  current  is  used  to  re¬ 
place  one  of  the  metrons.  The  calculation  of  fields  and 
currents  is  repeated  except  now  a  different  metron  is  re¬ 
placed  with  the  calculated  current.  This  is  repeated  un¬ 
til  sufficient  accuracy  is  obtained  or  there  cire  no  more 
metrons  to  update.  For  each  iteration,  it  was  found  that 
there  is  virtually  no  change  in  the  scattered  field  in  regions 
where  there  is  significant  energy.  The  iterative  scheme  is 
only  necessary  to  obtain  improved  accuracy  where  the 
fields  are  significantly  less  than  the  maximum.  There  are 
a  myriad  of  other  iterative  schemes  that  could  be  devel¬ 
oped;  the  one  presented  here  is  not  necessarily  the  opti¬ 
mum  one.  A  better  initial  selection  of  metrons  may  have 
eliminated  the  need  for  iteration.  However,  this  scheme 
does  illustrate  that  an  iterative  technique  can  be  used 
when  there  is  little  or  no  a  priori  knowledge  of  the  actual 
source  distribution. 

It  is  worth  mentioning  that,  in  principle,  an  iterative 
scheme  can  be  done  with  little  additional  computational 
overhead.  The  majority  of  CPU  time  is  spent  Ccdcu- 
lating  the  terms  G(r|r,)  for  pairs  of  surface  and  mesh 
points.  In  principle,  these  terms  just  need  to  be  calcu¬ 
lated  once  (since  they  depend  only  on  the  geometry  of 
the  grid)  and  then  multiplied  later  by  the  appropriate 
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metrons.  However,  in  practice  this  requires  the  storage 
of  2N^  terms,  which  represents  a  prohibitive  amount  of 
memory  for  Icirge  problems.  The  code  developed  for  this 
study  was  not  optimized,  and  the  Green’s  functions  were 
recalculated  for  each  iteration. 

To  calculate  fax-field  quantities  or  to  update  the 
metrons  in  the  iterative  scheme,  it  is  necessary  to  obtain 
the  surface  currents  from  the  MEI-derived  electric  field. 
The  currents  are  found  using 


-1  dE, 

jkoT]o  dn 


(23) 


where  n  is  normal  to  the  surface.  Given  the  mesh  struc¬ 
ture  and  the  fact  that  the  total  field  is  zero  at  the  surface, 
the  current  can  be  found  at  a  point  on  the  surface  using 


j  «  -1  E,{Av)  ^  -1 

jkorio  Av  jkoT]o  Av 

(24) 

where  E^^^(Av)  is  the  MEI-derived  scattered  field  at 
the  first  mesh  point  above  the  surface  (i.e.,  the  first  point 
above  the  surface  normal  and  on  the  surface  normal)  and 
E^^^iAv)  is  the  known  incident  field  at  the  same  point. 
However,  this  is  only  accurate  to  order  0(Av).  Since 
the  scattered  field  is  available  at  two  mesh  points  above 
the  surface,  the  surface  current  is  obtained  using  second- 
order  forward-differencing.  Thus,  the  normal  derivative 
is  found  using 

dE^  1 

^  ^  -  ^^^(2At;))  (25) 


where  E^  =  is  the  total  field. 

Our  goal  is  to  calculate  the  bistatic  radar  cross  section 
per  unit  length  for  a  plane  wave  incident  on  a  1-D  surface. 
This  is  found  using  [31] 


a{e„ei)  =  2it^  (26) 

liL 

where  p  is  the  distance  to  the  fax-field  observation  point, 
L  is  the  length  of  the  surface,  E  is  the  scattered  inten¬ 
sity,  li  is  the  incident  intensity,  Og  is  the  angle  of  obser¬ 
vation  me^lsured  from  the  vertical,  and  0,-  is  the  angle  of 
incidence  measured  from  the  vertical.  In  order  to  find 
the  bistatic  radar  cross  section,  it  is  necessary  to  convert 
the  surface  current  to  the  electric  field  in  the  far  field. 
Equation  (22)  could  be  used;  however,  given  the  distant 
location  of  the  observation  point,  the  large  argument  ap¬ 
proximation  of  the  Hankel  function  is  used.  Thus,  for  an 


incident  field  with  unit  magnitude,  we  obtain  the  cross 
section  using 

<T(es,ei)=  (27) 

ko  ^ 

■J  Js{ri)AEexp{jko{xi cos (9*  +  t/,-  sin E)) 

1  =  1 

V.  Results 

To  illustrate  the  results  of  the  method,  consider  a  sur¬ 
face  with  kgh  =  0.33  and  kol  =  2.83.  This  surface  has  an 
RMS  slope  angle  of  7  =  9.46°.  Fifty  surface  realizations 
of  length  L  =  80A  were  generated  and  enclosed  in  a  mesh 
with  Ax  =  An  =  A/16.  (This  mesh  spacing  is  typical  of 
that  reported  in  other  applications  of  the  MEI  method. 
It  was  found  that  using  a  finer  spacing  did  not  improve 
the  results,  while  a  coai-ser  mesh  yielded  slightly  worse 
results.)  Figure  6  shows  the  bistatic  radar  cross  section 
obtained  using  the  MEI  method  and  an  FDTD  technique 
[32]  that  is  known  to  be  accurate  for  this  type  of  sur¬ 
face.  Since  normal  incidence  was  used,  coherent  specular 
return  is  at  0°.  The  small  oscillations  in  the  curve  are 
due  to  the  use  of  a  finite  number  of  surface  realizations. 
Since  two  different  sets  of  surface  realizations  were  used, 
the  oscillations  do  not  coincide.  These  two  methods  show 
excellent  agreement  throughout  the  60  dB  range  of  the  re¬ 
sults.  The  MEI  method  results  were  obtained  after  three 
iterations  of  the  metrons. 

Figure  7  shows  the  cross  section  when  the  same  surface 
as  used  for  the  results  in  Fig.  6  is  illuminated  by  a  wave 
incident  at  45°.  All  aspects  of  the  calculation  (e.g.,  mesh 
structure,  taper  factor,  and  surface  realizations)  are  the 
same  as  before  except  five  iterations  were  used.  As  ex¬ 
pected,  the  coherent  specular  peak  is  at  45°.  The  two 
results  show  excellent  agreement  everywhere  except  to¬ 
ward  backward  grazing  angles. 

Figure  8  shows  the  results  for  a  surface  with  koh  =  1.72 
and  kol  =  7.31  that  is  illuminated  with  a  normally  inci¬ 
dent  wave.  Five  iterations  were  used.  Compared  to  the 
previous  surface,  the  correlation  length  has  increased  by 
a  factor  of  2^  and  the  RMS  surface  height  has  increased 
by  twice  this  same  factor.  This  change  in  parameters  ap¬ 
proximately  doubles  the  RMS  slope  angle  so  that  it  is 
now  18.43°.  Because  this  surface  is  fairly  rough,  it  pro¬ 
duces  no  clear  specular  return;  instead  incoherent  energy 
is  broadly  distributed  over  a  range  of  scattering  angles. 
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SCATTERING  ANGLE  [deg] 

Figure  6:  Radar  cross  section  as  a  function  of  scattered  angle  for  normal  illumination  of  a 
siudace  with  koh  =  0.33  and  kol  =  2.83. 


Again,  there  is  good  agreement  between  the  two  meth¬ 
ods  except  at  low  grazing  angles  where  the  MEI  method 
overpredicts  the  cross  section. 

The  three  results  presented  here  were  done  in  a  consis¬ 
tent  manner,  i.e.,  all  were  solved  using  the  same  approach. 
The  only  difference  between  runs  was  the  number  of 
iterations — no  other  modifications  were  attempted  to  try 
to  optimize  the  results  for  a  pcirticulax  geometry.  In  fact, 
there  is  every  reason  to  believe  that  the  results  presented 
here  are  not  optimum.  For  example,  a  better  metron 
for  use  in  the  initial  set  might  be  one  obtained  from  the 
physical  optics  equivalent,  ncimely,  Jj  «  2n  x  This 

could  be  used  instead  of  the  constant  metron  J}  given 
in  (21)  and  probably  would  be  appropriate  for  surfaces 
with  relatively  long  correlation  lengths.  Additionally,  by 
incorporating  the  x  and  y  components  of  the  wavenum¬ 
ber  rather  than  just  using  ko,  the  metrons  could  contain 
more  knowledge  of  the  incident  field.  There  may  also  be 
better  mesh  structures  than  the  one  presented  here. 


V.  Conclusions 

Although  the  limits  of  the  MEI  method  have  not  been 
completely  explored,  the  work  presented  in  this  paper  il¬ 
lustrates  that  the  method  has  the  potential  to  provide 
accurate  solutions  to  rough  surface  scattering  problems. 
The  problem  examined  in  this  study  required  a  large  dy¬ 
namic  remge.  Accuracy  was  achieved  using  an  iterative 
scheme.  Additional  accuracy  was  obtained  by  calculating 
the  surface  currents  using  second-order  forward  differenc¬ 
ing. 

Since  the  MEI  method  provides  a  solution  via  a  sparse 
and  relatively  small  matrix,  it  can  potentially  solve  rough 
surface  scattering  problems  in  three  dimensions.  In  that 
case  the  metrons  would  be  a  function  of  two  variables, 
rather  than  one,  and  it  may  be  much  harder  to  select 
“good”  metrons.  Therefore,  the  use  of  an  iterative  scheme 
in  three  dimensions  may  be  critical  for  successful  imple¬ 
mentation. 

The  TE  problem  has  not  been  discussed  but  is  solved  in 
a  manner  very  similar  to  that  presented  here.  The  solu¬ 
tion  for  penetrable  objects  involves  enclosing  the  material 
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Figure  7:  Radar  cross  section  as  a  function  of  scattered  angle  when  Oi  =  45°  for  a  surface 
with  fcfl/i  =  0.33  and  kol  =  2.83. 


interface  in  a  mesh — with  mesh  points  on  both  sides  of 
the  interface.  This  mesh  can  be  terminated  close  to  the 
interface,  but  the  solution  is  complicated  by  the  fact  that 
metrons  must  be  specified  in  terms  of  both  their  vedue  at 
the  surface  and  their  normal  derivative  at  the  surface. 

Several  other  problems  remciin  to  be  investigated.  For 
example,  if  a  surface  is  extremely  rough,  the  computa- 
tioncd  mesh  may  cross  itself.  In  other  words,  as  mesh 
points  are  placed  along  the  surface  normal,  some  points 
will  no  longer  have  monotonically  increasing  values  in  the 
X  direction  if  the  normal  directions  change  too  abruptly 
(such  as  in  a  narrow  valley).  This  type  of  highly  distorted 
mesh  does  not  permit  use  of  the  local  polar  approxima¬ 
tion  but  may  yield  to  analysis  by  using  the  MEI  method 
for  nodes  in  the  vicinity  of  the  mesh  overlap.  Other  re¬ 
maining  topics  include  determination  of  optimum  mesh 
structure  (both  in  terms  of  spacing  and  number  of  layers) 
and  best  selection  of  metrons. 
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Abstract 

WIREGRID  is  a  pre-processor  for  the  widely  used  pub¬ 
lic  domain  method  of  moments  computer  code  NEC2. 
WIREGRID  was  developed  for  modelling  metallic  struc¬ 
tures  using  a  mesh  of  wires.  In  particular,  it  has  pow¬ 
erful  facilities  for  controlling  the  density  of  the  mesh, 
permitting  the  same  basic  model  to  be  used  for  different 
frequencies,  with  the  code  automatically  generating  a 
suitable  mesh  for  a  particular  frequency.  WIREGRID  is 
not  limited  to  generating  wire  meshes:  it  can  also  be  use¬ 
fully  applied  to  true  thin-wire  modelling,  for  instance  for 
Yagi-Uda  and  log-periodic  antennas.  In  this  paper,  the 
basic  philosophy  of  the  code  is  described,  and  the  code 
capabilities  are  detailed.  Sufficient  detail  is  presented  so 
that  the  paper  can  serve  as  a  reference  for  code  users. 
The  meshing  algorithm  is  briefly  described.  An  example 
is  shown  of  the  automatic  mesh  generation  capabilities 
applied  to  an  automobile.  WIREGRID  is  available  in  the 
public  domain  via  anonymous  ftp;  details  of  the  ftp  sites 
and  hardware  requirements  are  given  in  the  appendices. 

1  Introduction 

NEC2  has  become  the  de-facto  standard  thin-wire 
method  of  moments  code  that  other  codes  are  evaluated 
against,  due  not  least  to: 

•  Its  powerful  facilities  for  general  thin-wire  mod¬ 
elling.  These  include  not  only  the  conventional  free- 
space  treatment  but  also  a  ground  treatment  (both 
an  approximate  reflection  coefficient  approach  and 
an  “exact”  —  in  a  numerical  sense  —  Sommerfeld 
treatment).  NEC2  can  also  exploit  symmetry,  and 
handle  non-radiating  networks  and  lumped  loads. 
This  brief  review  highlights  only  some  of  its  capa¬ 
bilities. 

•  The  extensive  documentation  available. 

•  The  widespread  availability  of  the  code. 


However,  anyone  who  has  used  NEC2  is  well  aware  that 
the  code,  powerful  as  it  is,  can  be  frustrating  and  time- 
consuming  to  use,  since  the  model  input  is  described  by 
a  list  of  coordinate  points.  It  is  extremely  easy  to  make 
an  error  in  such  an  input  file,  and  very  diflBcult  to  pick  it 
up  without  some  graphical  aid.  Furthermore,  when  such 
an  error  is  present,  the  code  will  often  run  and  produce 
results  that  are  plausible  but  incorrect.  Hence  a  number 
of  pre-processors  have  been  written  specifically  for  NEC2 
to  simplify  model  generation,  one  of  which  is  the  subject 
of  this  paper,  namely  WIREGRID,  which  we  are  mak¬ 
ing  available  in  the  public  domain.  We  assume  a  basic 
familiarity  with  NEC2  in  this  paper,  since  WIREGRID 
was  written  specifically  for  NEC2  model  generation. 

DIDEC  was  one  of  the  first  pre-processors  for  NEC2 
[1].  However,  IGUANA  (available  from  ACES  as  part  of 
the  NEEDS  package)  was  one  of  the  first  widely  available 
such  pre-processors;  it  has  been  generally  acknowledged 
as  a  useful  tool  but  with  a  number  of  deficiencies.  Work 
has  also  been  reported  by  other  researchers  on  graphical 
tools  for  preparing  NEC2  data  files;  this  journal  has  pub¬ 
lished  at  least  two  such  papers  [2,  3]  and  the  1993  ACES 
conference  saw  a  number  of  codes  of  varying  quality  and 
cost  being  demonstrated.  As  this  paper  went  to  print, 
another  South  Africa  group  working  at  the  University 
of  the  Witwatersrand  announced  a  commercial  software 
package  called  SIG  (Structure  Interpolation  and  Grid- 
ding),  which  offers  facilities  similar  to  or  exceeding  those 
in  WIREGRID  [4]. 

Our  pre-processor  WIREGRID  was  originally  devel¬ 
oped  for  a  specific  application  of  NEC2,  namely  mod¬ 
elling  metallic  structures  using  a  mesh  of  wires,  from 
whence  the  name.  In  particular,  it  has  very  powerful 
facilities  for  controlling  the  density  of  the  mesh,  permit¬ 
ting  the  same  basic  model  to  be  used  for  different  fre¬ 
quencies,  with  the  code  automatically  generating  a  suit¬ 
able  mesh  for  a  particular  frequency  (according  to  the 
A/10  rule,  for  instance).  Anyone  who  has  had  to  gener¬ 
ate  a  wire  grid  manually  or  even  using  a  pre-processor 
such  as  IGUANA  will  know  that  this  is  an  arduous  task: 
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specifying  the  end  point  coordinates  of  each  wire  seg¬ 
ment  of  the  mesh  (of  even  relatively  simple  structures) 
in  the  format  required  by  NEC2  is  an  error-prone  and 
exhausting  task  when  done  once;  and  if  total  number  of 
segments  has  to  be  changed  to  produce  a  finer  or  coarser 
discretization,  the  structure  must  almost  be  re-modelled 
from  scratch  for  each  change.  WIREGRID  waa  devel¬ 
oped  to  simplify  these  tasks  considerably.  The  initial 
number  of  coordinates  that  must  be  typed  in  by  hand 
is  typically  reduced  by  at  least  one  order  of  magnitude 
for  coarse  meshes,  and  two  or  even  three  orders  for  fine 
meshes.  Since  many  structures  can  be  modelled  by  a 
number  of  fiat  or  almost  flat  polygonal  surfaces,  it  was 
decided  that  a  pre-processor  was  needed  which  would 
take  as  input  just  the  vertices  of  these  polygons,  and 
which  would  then  generate  the  wiregrid  discretization 
automatically.  Changing  the  segment  density  then  only 
requires  a  few  keystrokes.  Once  the  user  is  satisfied  with 
the  model,  the  code  automatically  generates  the  NEC2 
input  file  with  the  geometry  cards  and  a  limited  number 
of  control  cards. 

WIREGRID  is  not  limited  to  generating  wire  meshes: 
it  can  be  usefully  applied  to  true  thin-wire  modelling, 
for  instance  for  Yagi-Uda  and  log-periodic  antennas,  and 
was  most  recently  used  for  a  project  involving  frequency 
selective  surface  modelling.  However,  most  of  its  features 
were  developed  with  wiregrid  modelling  in  mind,  since 
at  the  time  the  first  version  was  developed  (1990)  we 
were  involved  with  the  electromagnetic  characterization 
of  HE  and  VHP  whip  antennas  mounted  on  ground  vehi¬ 
cles  and  we  were  encountering  serious  problems  prepar¬ 
ing  the  large  data  sets  that  this  application  required. 
Closely  related  research  on  developing  a  parallel  version 
of  NEC2  was  also  in  progress  at  that  time  [5,  6];  due  to 
the  increased  problem  size  that  this  parallel  version  of 
NEC2  could  handle,  WIREGRID  also  proved  very  useful 
with  this  work. 

The  code  was  originally  developed  by  the  first  author 
in  close  liaison  with  the  second  author,  and  has  sub¬ 
sequently  been  extensively  used  for  both  vehicle  EMC 
modelling  under  contract  (as  outlined  above)  and  teach¬ 
ing  (at  undergraduate  and  postgraduate  level,  as  well  as 
for  short  courses  on  NEC2).  The  idea  underlying  the 
code  is  similar  to  that  reported  by  Najm  [3],  but  the  im¬ 
plementation  is  rather  different.  WIREGRID  is  highly 
graphically  orientated,  whereas  Najm’s  code  is  not.  It 
also  automates  a  number  of  functions  that  Najm’s  code 
requires  manual  intervention  to  implement;  see  for  ex¬ 
ample  the  discussion  regarding  the  motor  vehicle  [3,  Fig. 
5].  Despite  the  conceptual  similarities  of  the  two  codes, 
WIREGRID  was  developed  independently  and  at  about 
the  same  time  [7]. 

As  with  most  complex  computational  electromagnetic 
codes,  while  a  pre-processor  such  as  WIREGRID  greatly 
simplifies  the  task  of  the  user,  a  knowledge  of  the  input 


requirements,  and  the  basic  functioning,  of  the  code  is 
still  highly  advisable.  The  user  may  need  to  append  ad¬ 
ditional  control  statements  to  the  input  file  and  may  even 
want  to  change  some  of  the  geometry  data  afterwards. 

2  Using  the  program 

2.1  Structure  Modelling 

The  tasks  required  by  the  user  are  the  following:  The 
structure  must  first  be  divided  into  fairly  flat  polygonal 
surfaces,  which  are  called  major  elements  (see  Fig.  1). 
WIREGRID  can  handle  up  to  160  such  major  elements. 
Each  major  element  does  not  need  to  lie  completely  in 
a  2D  plane,  but  severely  “twisted”  elements  are  really 
beyond  the  capability  of  the  code.  The  frequent  visual 
feedback  available  from  the  program  will  permit  the  user 
to  identify  unsuitable  elements  and  rectify  them.  As 
mentioned  in  the  Introduction,  the  code  was  developed 
for  structures  consisting  mainly  of  flat  surfaces;  it  is  not 
really  suited  modelling  curved  surfaces,  such  as  an  air¬ 
craft  fuselage.  The  coordinates  of  each  major  element’s 
vertices,  called  major  nodes,  must  be  determined  next 
(see  Fig.  1).  The  simplest  major  element  is  defined  by 
two  major  nodes,  which  would  model  a  single  straight 
wire.  The  most  complex  major  element  may  be  defined 
by  a  maximum  of  10  major  nodes. 

When  generating  a  mesh,  it  is  important  to  ensure 
that  meshes  in  adjacent  elements  are  properly  electri¬ 
cally  connected  across  the  boundary.  (NEG2  users  will 
be  aware  of  a  particularly  insidious  fault  caused  by 
two  wires  whose  end-points  that  are  closely  located  and 
should  be  physically  connected,  but  are  treated  by  the 
code  as  unconnected.)  We  will  use  the  term  conformal¬ 
ity  in  this  paper  to  describe  a  properly  connected  mesh. 
WIREGRID  will  ensure  conformality  provided  that  ma¬ 
jor  nodes  of  adjacent  major  elements  coincide  on  the 
boundary  that  the  elements  share.  These  requirements 
mean  that  when  vertices  on  the  boundary  between  two 
adjacent  polygons  are  not  coincident,  extra  major  nodes 
for  one  or  both  of  them  must  be  defined.  These  extra 
nodes  must  coincide  with  the  unmatched  vertices  (see 
Fig.  2).  The  same  technique  can  be  used  for  ensuring 
that  a  whip  antenna  is  properly  grounded  to  the  vehi¬ 
cle  frame.  Superfluous  nodes  (for  example,  a  node  in 
the  middle  of  an  element  side  when  there  are  no  confor¬ 
mality  problems)  do  not  cause  any  problem  within  the 
program  (providing  of  course  that  all  elements  sharing 
that  side  contain  the  extra  node).  However,  the  wire 
mesh  generated  may  be  distorted  by  their  presence,  and 
it  is  recommended  that  such  unnecessary  major  nodes 
be  avoided. 

An  additional  geometric  requirement  for  the  major  el¬ 
ements  is  that  they  must  be  convex,  i.e.  a  straight  line 
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major  nodes 


Figure  1:  Arbitrary  structures  can  be  built  up  with 
flat  polygonal  surfaces.  The  example  shows  a  simplified 
model  of  a  generic  “hatchback”  motor  car. 


joining  any  two  points  in  the  polygon  must  still  be  en¬ 
closed  by  the  polygon  boundaries.  The  reason  for  this 
requirement  is  that  it  greatly  simplifies  the  task  of  the 
mesh  generator  in  WIREGRID.  It  can  be  adhered  to 
easily,  since  any  non-convex  polygon  can  be  divided  into 
two  or  more  convex  parts. 

A  number  tag  starting  from  0  is  assigned  by  WIRE- 
GRID  to  each  major  element.  It  is  often  helpful  to  keep 
a  rough  sketch  of  the  structure,  indicating  all  the  major 
elements  with  their  number  tags.  In  the  NEC2  input 
file  that  WIREGRID  generates,  all  wires  belonging  to  a 
certain  major  element  will  be  tagged  with  this  number. 
When  certain  wires  are  to  be  loaded,  excited  or  even 
altered  by  hand,  this  will  help  to  locate  them. 

2.2  Discretization 

The  nominal  length  in  of  the  wire  segments  in  the  mesh 
may  be  changed  at  any  time  during  the  structure  mod¬ 
elling  process.  The  new  discretization  is  performed  im¬ 
mediately.  The  nominal  segment  length  is  a  key  concept 
in  the  code:  each  straight  wire  is  divided  into  p  segments, 
defined  by 

p  =  intf^ -I- 0.99999V  (1) 


where  L  is  the  length  of  the  wire  section  and  int{a:)  is 
the  integer  operator,  taking  the  integer  part  of  x. 

All  straight  wires  belonging  to  a  major  element  defined 
by  more  than  two  major  nodes  are  divided  into  segments 
of  equal  length.  Following  the  “same  surface  rule”  [8], 
these  wires  are  modelled  with  a  radius  of 

When  the  wire  is  a  major  element  defined  by  two 
nodes,  the  “same  surface  rule”  is  irrelevant,  and  the  user 
must  specify  radii  for  the  two  end  points,  which  is  in¬ 
dependent  of  the  discretization.  (The  code  still  offers  a 
radius  of  as  the  default,  although  this  will  of  course 
almost  certainly  be  incorrect  in  this  case).  The  number 
of  segment  divisions  on  the  wire  is  still  defined  by  (1).  In 
addition  to  simple  wires  of  uniform  radius  (correspond¬ 
ing  to  a  GW  NEC2  card),  WIREGRID  can  also  handle 
tapered  wires,  generating  the  necessary  GW  and  GC  cards. 


2.3  A  brief  overview  of  the  meshing  al¬ 
gorithm 


Figure  2:  To  ensure  that  the  mesh  within  each  major  el¬ 
ement  is  electrically  continuous  across  the  common  side, 
extra  major  nodes  must  sometimes  be  added  as  shown 
here. 


A  detailed  description  of  the  meshing  algorithm  is  be¬ 
yond  the  scope  of  this  paper.  The  following  description 

^One  reviewer  expressed  surprise  at  this  radius,  but  Ludwig’s 
paper  does  indeed  verify  this  rule.  This  rule  actually  refers  to  each 
polarization,  so  as  Ludwig  comments,  “twice  surface  area  rule” 
might  be  a  better  name.  Trueman  and  Kubina  give  the  same 
formula  for  wire  radius  [2,  p.21].  Wires  of  this  radius  are  on  the 
verge  of  being  too  thick  for  the  thin-wire  treatment  in  NEC2  (the 
extended  kernel  option  does  not  help  with  a  mesh,  since  it  cannot 
be  used  for  segments  with  a  complex  junction),  but  experience  has 
produced  generally  acceptable  results  following  this  rule. 
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is  intended  to  give  only  a  flavour  of  the  meshing  algo¬ 
rithm.  Figure  3  shows  the  significant  features  of  the 
meshing  process,  and  will  be  referred  to  during  this  dis¬ 
cussion.  Five  major  nodes  are  shown  by  •;  only  four  are 
actually  required  for  this  major  element,  but  the  fifth 
would  be  required  if  another  major  element  were  to  have 
a  major  node  at  this  point.  (See  for  instance  Figure  2 
for  such  an  example.)  At  the  risk  of  tedium,  we  reiterate 
that  the  coordinates  of  these  major  nodes  are  all  that  the 
user  need  enter. 

•  The  sides  of  the  polygonal  are  divided  up,  aa  close 
to  the  specified  nominal  segment  length  aa  possible. 
These  generate  minor  nodes  along  the  sides.  The 
minor  nodes  are  represented  by  o  in  Figure  3.  These 
minor  nodes  are  then  used  to  the  basis  for  the  mesh 
as  follows: 

•  The  vertex  with  included  angle  closest  to  90°  is  then 
found. 

•  The  sides  of  the  polygon  are  then  numbered  counter¬ 
clockwise.  Assume  that  the  polygon  is  close  to  rect¬ 
angular;  this  will  then  generate  sides  1  to  4. 

•  Starting  at  the  first  point  on  side  1  where  a  minor 
node  should  be  located,  the  code  constructs  vertical 
cross-wires  from  side  1  to  side  3,  endeavouring  to 
keep  the  cross-wires  as  close  as  possible  to  parallel 
to  each  other.  (These  wires  are  coloured  red  on  the 
screen). 


Side  3 


Side  4 


Side  1 


Figure  3:  A  rectangular  major  element  illustrating  the 
meshing  algorithm  outlined  in  Section  2.3.  See  text  for 
discussion. 
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[3]  are  not  good  meshes  for  NEC2:  with  a  significant 
wire  radius  they  are  very  likely  to  generate  match 
points  inside  another  wire’s  volume;  the  NEC2  user 
guide  specifically  cautions  against  this  [9,  p.3-6]). 
Users  should  try  to  avoid  entering  major  elements 
with  such  geometries.  Automatically  detecting  such 
violations  of  the  NEC2  rules  and  guidelines  is  a  fea¬ 
ture  we  would  like  to  add  but  have  not  yet;  we  will 
comment  on  this  later. 


•  Once  these  cross-wires  have  been  constructed,  the 
cross-wires  connecting  sides  2  and  4  are  constructed 
in  a  similar  fashion.  (These  are  coloured  blue). 
The  mesh  thus  created  is  shown  in  Figure  3.  The 
mesh  is  represented  by  the  dashed  lines  in  the  figure. 
In  general,  the  horizontal  cross-wires  are  individual 
wires  running  between  each  vertical  cross- wire,  with 
potentially  slightly  different  lengths  and  directions. 
However,  due  to  the  geometrical  simplicity  of  this 
example,  all  cross-wires  in  Figure  3  are  uniform. 

•  The  code  has  a  number  of  features  to  try  to  gener¬ 
ate  a  “good”  mesh  for  irregularly  shaped  polygons; 
for  instance,  a  point  on  a  side  may  sometimes  be 
skipped  if  a  “better”  mesh  can  be  constructed  by 
using  the  next  point  on  a  side  as  the  end-point  for  a 
cross- wire.  (“Good”  and  “better”  refer  to  the  angles 
that  the  cross- wires  make  with  the  element  sides.) 
The  intersections  of  the  “blue”  and  “red”  wires  may 
also  be  perturbed  slightly  with  the  aim  of  making 
the  included  angles  between  the  wires  as  close  as 
possible  to  90°.  (Neither  of  these  mesh  improve¬ 
ment  stages  are  required  in  the  mesh  of  Figure  3.) 

•  The  code  tries  to  avoid  creating  long,  thin  triangular 
meshes.  (The  triangular  examples  shown  by  Najm 


2.4  Generation  of  the  NEC2  input  file 

WIREGRID  supports  two  data  formats  for  the  NEC2 
input  file: 

•  Fixed  length  data  fields.  All  numbers  are  separated 
by  spaces,  where  all  integers  occupy  fields  three  to 
five  characters  long,  and  floating  point  numbers  oc¬ 
cupy  10  character  long  fields.  This  is  the  format 
documented  in  the  user  guide  [9]. 

•  Variable  length  data  fields.  All  numbers  are  sepa¬ 
rated  by  commas.  A  number  of  newer  versions  of 
NEC2  also  support  this  format. 

When  satisfied  with  the  model,  the  first  part  of  the  NEC2 
input  file  containing  all  the  geometric  data  may  then  be 
generated  on  command. 

WIREGRID  also  provides  for  the  generation  of  control 
statements  which  allow  maximum  coupling  computation 
between  voltage  sources,  and  radiation  pattern  genera¬ 
tion  at  specified  frequencies.  Only  wire  segments  be¬ 
longing  to  major  elements  defined  by  two  nodes  may  be 
excited  for  this  purpose.  (Of  course,  the  experienced  user 
can  edit  the  NEC2  data  file  created  by  WIREGRID  to 
define  more  complex  antennas).  Furthermore,  the  radi¬ 
ation  pattern  cards  generated  are  only  for  cross  sections 
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defined  by  angles  of  constant  azimuth  or  elevation  angles 
of  the  full  three  dimensional  patterns.  Since  it  is  fairly 
ea.sy  to  enter  the  control  statements  into  the  NEC2  input 
file  by  hand,  only  these  very  common  type  of  requests 
are  catered  for  by  WIREGRID. 

Any  additional  command  lines  or  alterations  to  the 
geometry  data  may  be  added  afterwards,  with  the  user’s 
editor  of  choice.  For  this  purpose,  the  following  infor¬ 
mation  on  the  structure  of  the  geometry  data  section  is 
important; 

•  Each  tag  number  in  the  file  refers  to  the  number 
of  the  major  element  wherein  the  wire  segment  is 
located.  This  means  that  wire  segments  belonging 
to  major  elements  defined  by  more  than  two  major 
nodes  do  not  have  unique  tag  numbers,  and  can¬ 
not  be  referred  to  by  the  control  statements  at  the 
end  of  the  file,  unless  the  tag  number  is  changed 
by  the  user  afterwards.  This  is  also  the  reason  why 
WIREGRID  only  allows  excitation  of  segments  in 
two-node  major  elements  in  its  control  statement 
generation  options. 

•  Geometry  data  for  the  polygonal  boundaries 
(yellow^)  between  the  major  elements  are  entered 
first  into  the  file,  starting  with  the  boundaries  of 
major  element  number  0,  and  then  the  two  node 
elements  are  added. 

•  This  is  followed  by  data  statements  describing  the 
internal  meshes  (where  relevant).  These  statements 
are  divided  into  sets  for  each  polygonal  major  el¬ 
ement,  starting  with  the  element  with  the  lowest 
tag  number.  In  addition,  the  internal  mesh  data 
statements  for  each  polygonal  major  element  are 
subdivided  into  a  subset  of  statements  describing 
long,  straight,  segmented  wire  sections  (red),  and  a 
subset  describing  short,  single  segment  wire  sections 
(blue). 


2.5  The  user  interface 

Fig.  4  shows  the  main  menu  screen  of  WIREGRID.  An 
option  is  selected  by  moving  the  cursor  to  the  required 
one  and  hitting  return.  Should  the  wrong  option  be 
selected  accidentally,  hitting  escape  at  any  time  while 
WIREGRID  is  prompting  for  input  information  returns 
the  user  to  the  main  menu,  without  any  action  being 
taken^. 

^The  colours  mentioned  are  used  in  the  graphical  display  of  the 
model  (see  Section  2.3). 

^Editing  major  element  coordinates,  or  selecting  wire  segments 
for  excitation,  are  exceptions  to  this  rule.  In  these  cases  hitting 
escape  ends  the  editing  or  selecting  session.  An  element  thus  in¬ 
correctly  created  is  removed  by  deleting  the  major  element. 


WIRE  GRID  MODELLING  FOR  NEC 


Load  a  model . 

Specify  nominal  segment  length. 
Add  a  major  element. 

Edit  a  major  element. 

Delete  a  major  element. 

Edit  a  maj  or  node . 

Translate  major  elements. 

Rotate  major  elements. 

Mirror  major  elements. 

Check  number  of  segments . 

View  a  major  element. 

View  model . 

Save  model . 

Create  input  file  for  NEC. 

Quit . 


Use  <t>  and  <J.>  to  choose  an  option 
Press  <ENTER>  to  execute  the  option 


Figure  4:  WIREGRID’s  main  menu  screen 


2.5.1  Creating  a  new  model 

When  a  new  model  is  started,  the  nominal  segment  length 
must  be  specified  before  any  major  elements  can  be  en¬ 
tered.  Experience  has  taught  many  new  users  have  prob¬ 
lems  with  the  code  at  this  initial  stage,  not  realizing  that 
a  nominal  segment  length  is  required.  Any  reasonable 
value  suffices;  as  already  discussed  the  nominal  segment 
length  can  be  changes  at  any  time  during  the  modelling 
process.  The  first  major  element  can  then  be  entered.  If 
a  model  is  already  loaded,  and  the  user  wants  to  start 
a  new  model,  it  is  easiest  to  quit  and  run  the  program 
again. 


2.5.2  Deleting  major  elements 

A  major  element  is  deleted  by  selecting  the  “Delete  a 
major  element”  option.  It  is  important  to  note  that  all 
elements  with  tag  numbers  higher  than  the  one  deleted 
will  be  assigned  a  new  number,  one  less  than  before. 
When  more  than  one  element  must  be  deleted,  it  is  there¬ 
fore  best  to  start  with  the  highest-numbered  element  and 
work  downwards.  This  eliminates  the  need  to  continu¬ 
ously  recalculate  the  tag  number  of  the  elements  still  to 
be  deleted. 
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2.5.3  Adding,  editing  and  saving  models 

An  existing  model  may  be  loaded  from  disk,  together 
with  the  last  discretization  used.  The  user  may  there¬ 
fore  start  to  edit  major  elements  or  add  new  ones  imme¬ 
diately  after  loading. 

When  a  new  major  element  is  defined,  WIREGRID 
prompts  for  the  element’s  number  tag  and  the  number  of 
nodes  defining  the  element.  The  default  value  of  the  lat¬ 
ter  is  two,  and  replacing  it  by  larger  numbers  (up  to  10), 
will  result  in  more  complex  elements.  WIREGRID  sub¬ 
sequently  presents  the  user  with  a  table  of  major  node 
coordinates,  with  default  values  corresponding  to  the  re¬ 
spective  nodes  of  the  previous  major  element^.  The  cur¬ 
sor  is  moved  around  in  the  table  with  the  left  and  right 
tab  keys,  and  with  the  up  and  down  arrow  keys,  to  edit 
the  various  coordinates.  Coordinate  editing  is  followed 
by  specification  of  the  end  point  radii  for  two-node  ma¬ 
jor  elements.  (The  default  radius  of  ^n/27r  has  already 
been  discussed  in  Section  2.2). 

The  default  values  for  the  latter  is  ^„/27r,  the  value 
used  for  the  polygonal  major  element  wire  sections  at 
that  stage.  (This  is  of  course  probably  wrong,  since  two- 
node  element  frequently  represent  wire  antennas,  but  it 
is  trivial  to  enter  the  correct  radius). 

Editing  an  existing  major  element  is  the  same  as  defin¬ 
ing  a  new  one,  except  that  the  default  values  for  the 
number  of  nodes,  coordinates  and  end  point  radii,  are 
the  actual  values  specified  previously. 

Models  can  be  combined  by  loading  an  existing  model 
from  disc  at  any  time  and  merging  it  with  the  current  one 
in  memory,  helping  the  user  to  model  the  structure  in  an 
organized  and  structured  fashion.  An  example  would  be 
the  hull  and  turret  of  a  ground  vehicle,  which  are  conve¬ 
niently  modelled  separately  and  then  combined.  When 
the  models  are  merged,  WIREGRID  ensures  that  ele¬ 
ments  from  the  different  models  with  common  edges  are 
correctly  picked  up,  ensuring  that  the  mesh  is  properly 
conformal. 

A  range  of  subsequent  major  elements  may  be  spa¬ 
tially  translated  by  adding  a  vector  to  all  coordinates 
involved.  Alternatively,  they  can  be  rotated  around  an 
axis  (specified  by  the  user  as  a  line  through  two  three- 
dimensional  Cartesian  points)  by  a  certain  angle,  or  mir¬ 
rored  in  a  plane  (specified  by  the  user  as  the  spherical 
coordinate  vector  from  the  origin  to  the  nearest  point  on 
the  plane  surface).  In  the  process  the  elements  can  be  ef¬ 
fectively  moved  to  these  new  positions,  or  new  elements 
can  be  created  in  this  way.  This  helps  the  user  to  build 
up  structures  with  symmetrical  properties  fairly  quickly. 
Note  that  these  “Translate” ,  “Mirror”  and  “Rotate”  op¬ 
tions  do  not  generate  associated  structure  geometry  lines 
(GM,  GX,  GR-  “cards”  respectively)  [9,  p.  16]  in  theNEC2 

■‘This  eliminates  to  a  certain  extent  the  need  to  retype  node 
coordinates  common  to  the  previous  element. 


input  file,  but  physically  generate  new  or  altered  major 
elements. 

The  model  can  be  saved  at  any  stage  by  specifying  an 
appropriate  file  name.  It  is  recommended  that  this  be 
done  often,  since  a  major  element  may  be  accidentally 
deleted,  or  an  unforeseen  run-time  error  may  occur.  The 
latter  is  unfortunately  not  infrequent,  since  it  is  very 
difficult  to  foresee  all  errors  that  may  occur,  in  particular 
with  incorrect  user  input.  Provided  that  the  model  has 
been  recently  saved,  it  is  simple  to  reload  WIREGRID, 
load  the  model  and  then  continue  building  the  model. 

2.5.4  Entering  numbers 

When  typing  in  a  number,  or  editing  one,  the  Backspace, 
Delete  and  the  Left  and  Right  arrow  keys  can  be  used 
in  the  usual  manner.  The  default  editing  mode  is  always 
the  “Insert”  mode,  but  it  can  be  toggled  to  the  “Over¬ 
write”  mode  with  the  Insert  key.  Ctrl  +  Left  or  Right 
arrows  will  move  the  cursor  to  the  end  or  beginning  of 
the  number  being  edited. 

2.5.5  Inspecting  the  model 

The  entire  model  can  be  inspected  at  any  time  by  se¬ 
lecting  the  “View  modef’-option,  or  a  single  major  ele¬ 
ment  may  be  inspected  by  selecting  the  “View  a  major 
element”-option  instead. 

The  boundaries  of  the  major  elements  are  displayed 
in  yellow,  while  wires  generated  by  the  discretizer  are 
displayed  in  red  and  blue.  These  have  been  discussed  in 
§2.3. 

The  view  of  the  model  can  be  changed  by  zooming 
in  or  out  with  the  +  and  —  keys,  by  rotating  it  with 
the  arrow  keys,  and  by  translating  it  sideways  or  up  and 
down  with  the  L,  R,  U,  and  D  keys. 

Any  major  element  can  be  highlighted  by  hitting  S  and 
specifying  the  major  element  number.  When  viewing 
a  single  major  element,  this  can  be  used  to  add  other 
major  elements  to  the  picture.  This  feature  is  the  only 
one  that  does  not  work  with  a  monochrome  graphics 
adapter,  since  the  highlighted  (white)  colour  cannot  be 
distinguished  from  the  other  colours. 

2.5.6  Excitation  of  wire  segments 

In  the  final  generation  of  the  NEC2  input  file,  the  user 
may  specify  the  excitation  of  some  wire  segments  for 
maximum  coupling  computation  and  radiation  pattern 
generation  (see  Section  2.4). 

For  maximum  coupling  calculations,  up  to  five  seg¬ 
ments  may  be  excited,  and  for  an  antenna  array  pattern 
calculation,  up  to  ten  segments  may  be  excited. 
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Figure  5;  A  passenger  motor  car  showing  only  the  major 
elements.  There  are  41  segments  in  this  model.  The 
model  was  generated  using  symmetry,  hence  the  wires 
on  the  centreline. 


A  list  of  selected  segments  (specified  by  their  major  el¬ 
ement  tag  number,  and  the  segment  number®)  is  shown 
near  the  top  of  the  computer  screen.  Below  that,  infor¬ 
mation  on  all  two-node  major  elements  is  displayed  for 
one  element  at  a  time,  while  the  user  is  “paging  through” 
the  list  with  the  left  and  right  arrow  keys.  In  order  to 
select  a  wire  segment,  the  user  “pages”  to  the  particular 
two-node  major  element,  and  increments  or  decrements 
the  segment  number  with  the  up  and  down  arrow  keys. 
The  insert  key  is  finally  used  to  select  a  segment  for 
excitation,  whereupon  WIREGRID  will  prompt  for  the 
excitation  voltage,  and  the  segment  will  be  added  to  the 
list  of  excited  segments.  The  delete  key  is  used  to  remove 
segments  from  this  list. 

We  have  found  from  experience  that  users  often  have 
trouble  using  this  particular  feature  of  WIREGRID, 
partly  because  their  is  no  visual  feedback  to  confirm  the 
location  of  the  sources.  Since  it  is  relatively  easy  to 
manually  append  the  necessary  source  commands  to  the 
NEC2  input  deck  created  by  WIREGRID,  many  users 
elect  to  do  this. 


2.6  An  example 

In  Figures  5  and  6  we  show  the  WIREGRID  mesh  for 
a  motor  vehicle,  very  similar,  although  not  identical,  to 
the  generic  hatchback  shown  in  Figure  1. 


®  Segment  numbers  start  at  1,  denoting  the  segment  joined  to 
node  0. 
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Figure  6:  The  same  passenger  motor  car  with  a  nominal 
segment  length  of  0.1m.  (Using  the  A/10  rule,  this  model 
would  be  valid  up  to  around  300  MHz).  The  number  of 
segments  is  now  3008.  It  must  be  emphasized  that  only 
one  key  stroke  was  required  to  generate  this  mesh  from 
the  model  shown  in  Figure  5.  Note  that  the  code  does 
not  do  hidden  line  removal.  (A  very  powerful  computer 
would  be  required  to  run  this  model:  the  discretization 
shown  here  is  to  emphasize  the  capabilities  of  WIRE- 
GRID). 

3  Conclusions 

We  have  found  that  WIREGRID  has  greatly  boosted  our 
productivity  for  generating  NEG2  input  decks  for  simu¬ 
lations  involving  wire  mesh  modelling.  WIREGRID  has 
a  number  of  limitations,  many  of  which  we  have  outlined 
in  this  paper,  but  there  are  some  other  fundamental  lim¬ 
itations  that  potential  users  should  be  aware  of.  WIRE- 
GRID  does  not  check  the  validity  of  the  structures  that 
it  generates  with  respect  to  all  the  NEG2  model  gener¬ 
ation  rules  —  admirably  summarized  by  Trueman  and 
Kubina  [2] .  Particularly  when  the  structure  involves  thin 
triangular  elements,  it  is  quite  possible  that  the  WIRE- 
GRID  mesh  may  violate  some  of  the  validity  conditions. 
Of  course,  its  output  could  easily  be  run  through  such  a 
validity  checker,  but  at  the  time  of  writing,  Trueman  and 
Kubina’s  code,  GHEGK,  is  still  not  generally  available 
due  to  export  restriction  problems.  We  would  like  to  ex¬ 
tend  WIREGRID  to  incorporate  at  least  some  of  these 
validity  checks.  We  would  also  like  to  add  a  facility  to 
read  in  and  display  NEC2  data  decks  in  WIREGRID;  it 
does  not  presently  have  such  a  facility.  However,  funding 
for  this  work  is  very  restricted  at  present  and  neither  fa¬ 
cility  is  likely  to  implemented  in  the  near  future.  Finally, 
while  we  have  found  the  user  interface  to  be  very  func¬ 
tional,  it  does  not  have  the  latest  Windows-style  icons, 
dialogue  boxes,  hot-keys  etc. 
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Appendices 

A  Obtaining  and  Installing  the 
Code 

The  code  is  available  on  the  NEC  Archives  ftp  site.  The 
ftp  site  name  is  ftp.netcom.com,  and  the  code  is  in 


pub/rander/NEC  directory  as  wiregrid.zip.  The  code 
was  compressed  using  the  widely  used  DOS  shareware 
utility  PKZIP;  use  PKUNZIP  to  decompress.  Details  on 
downloading  codes  from  anonymous  ftp  sites  may  be 
found  in  any  of  the  many  guides  to  the  Internet  now 
appearing.  To  summarize,  proceed  as  follows: 

•  Enter  ftp  ftp.netcom.com 

•  The  system  will  respond  with  Name:  (or  similar 
message);  reply  with  anonymous. 

•  The  system  will  then  prompt  for  a  password;  it  is 
customary  to  enter  your  Internet  (e-mail)  address. 

•  Enter  cd  /pub/rander/NEC. 

•  then  binary  (this  is  actually  the  default  for  this 
site). 

•  and  finally  get  wiregrid.zip. 

•  To  exit  the  ftp  session,  quit. 

There  may  be  an  appreciable  delay  (possibly  a  number 
of  seconds)  between  entering  a  command  and  obtaining 
a  response,  especially  if  you  are  accessing  the  site  from 
outside  the  USA.  South  African  users  should  note  that 
wiregrid.zip  is  mirrored  on  our  local  Stellenbosch  ftp 
site,  viz.  ftp.sun.ac.za,  and  should  download  from 
this  site  in  preference  to  ftp.netcom.com. 

If  PKUNZIP  cannot  unzip  the  file,  check  that  you  have 
a  sufficiently  recent  version  (the  index  entry  for  WIRE- 
GRID  details  the  version  used);  if  this  is  not  the  prob¬ 
lem,  you  may  have  done  a  ftp  transfer  in  ASCII  mode  en 
route  to  your  PC.  Your  local  Internet  guru  will  be  able 
to  advise! 

For  installation,  all  files  on  the  distribution  disk  should 
be  copied  to  a  directory  of  the  user’s  choice.  When 
wiregrid.exe  is  run,  it  will  prompt  for  the  directories 
where  the  wire  grid  models  and  the  NEC2  input  files 
are  to  be  stored  respectively.  Several  example  files  are 
included  with  the  zip  file.  The  user  may  change  these 
directories  during  a  session,  and  these  changes  then  be¬ 
come  the  defaults  for  the  next  session.  Whenever  a  direc¬ 
tory  which  does  not  exist  is  specified,  it  will  be  created 
on  request. 

B  Hardware  requirements 

WIREGRID.EXE  runs  under  DOS  on  an  IBM  com¬ 
patible  PC.  The  program  is  written  in  TurboPascal, 
and  makes  extensive  use  of  the  TurboPascal  graphics 
calls.  As  a  result,  it  will  unfortunately  not  be  a  sim¬ 
ple  matter  to  port  to  another  hardware  platform  using 
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a  different  operating  system  (such  as  an  Apple  Macin¬ 
tosh).  Although  a  8087/80287/80387  mathematical  co¬ 
processor  would  be  advantageous,  it  is  not  a  prerequisite. 
(The  code  automatically  detects  the  presence  of  a  co¬ 
processor).  A  colour  graphics  adapter  and  a  colour  dis¬ 
play  monitor  (or  a  monochrome  display  which  can  differ¬ 
entiate  between  colours  using  various  shades  of  grey)  are 
required  to  fully  exploit  the  graphical  features  of  WIRE- 
GRID,  but  almost  all  of  these  will  still  work  satisfacto¬ 
rily  with  a  monochrome  graphics  adapter.  WIREGRID 
supports  the  GGA,  EGA,  IBM8514,  Hercules,  ATT400, 
VGA  and  the  PC3270  graphics  adapters. 

C  Recent  Improvements 

This  paper  describes  WIREGRID  Version  3.  The  latest 
version,  4.00,  incorporates  some  useful  new  features  for 
viewing  the  model.  These  are  : 

•  The  model  can  now  be  plotted  as  if  viewed  along  on 
of  the  three  axes,  X,  Y  or  Z.  (Actually,  there  are 
six,  not  three,  options  here,  since  x  and  X  rotate 
the  view  by  180°.) 

•  The  user  can  now  toggle  perspective  on  and  off.  The 
default  is  to  show  the  model  in  perspective,  but 
there  are  cases  where  it  is  more  convenient  to  re¬ 
move  this  facility,  such  as  when  the  model  is  viewed 
along  an  axis. 

•  The  display  can  now  be  switched  to  monochrome 
mode  on  colour  screens.  (There  is  no  effect  on 
monochrome  screens).  This  is  useful  for  doing 
screen-prints  and  screen-captures;  for  example,  the 
screen-capture  utility  that  was  used  to  generate  Fig¬ 
ures  5  and  6  in  this  paper  produces  satisfactory  im¬ 
ages  from  a  monochrome  display,  but  not  from  a 
coloured  one. 

D  Further  Information 

Users  are  welcome  to  contact  the  second  author  at  the 
address  given  on  this  paper  to  report  bugs  or  make  sug¬ 
gestions,  but  as  with  most  public  domain  software,  we 
unfortunately  cannot  offer  extensive  user  support.  This 
paper  serves  as  the  user  manual  at  present;  no  further 
documentation  is  available. 

The  standard  software  disclaimer  applies:  no  warran¬ 
tee  is  given  or  implied  that  all  errors  have  been  elimi¬ 
nated  from  the  program  WIREGRID,  and  the  authors 
and  the  University  of  Stellenbosch  accept  no  responsibil¬ 
ity  for  any  losses,  damages  etc.  resulting  from  the  use  of 
the  program. 


The  authors 

Cornelis  (Nelis)  F.  du  Toit  was  born  in  South  Africa  on  May 
3,  1962.  He  received  the  B.Eng.,  M.Eng.  and  Ph.D.  from 
the  University  of  Stellenbosch  in  1984,  1986  and  1992  respec¬ 
tively.  His  doctoral  dissertation  was  entitled  “The  computa¬ 
tion  of  electromagnetic  fields  from  cylindrical  near-field  mea¬ 
surements  at  non-asymptotic  distances  from  an  antenna”. 
He  was  a  Research  Assistant  at  the  University  from  1988- 
1992,  and  worked  at  Plessey-Tellumat  in  Cape  Town  during 
1993  as  a  contractor.  Nelis  and  his  family  moved  to  New 
Zealand  in  January  1994,  and  he  is  presently  working  for 
DelTec  in  Wellington.  The  work  reported  in  this  paper  was 
performed  whilst  Nelis  was  at  the  University  of  Stellenbosch. 
His  research  interests  include  electromagnetic  field  theory, 
antenna  design,  near-field  measurement  techniques  and  nu¬ 
merical  methods  in  antenna  engineering.  Nelis  is  a  keen  am¬ 
ateur  astronomer  and  an  accomplished  pianist;  he  also  enjoys 
listening  to  music.  Nelis  is  married  to  Elize  (a  musician)  and 
they  have  one  child. 

David  B.  Davidson  (IEEE  Member)  was  born  in  London, 
England  on  April  13,  1961.  He  grew  up  in  South  Africa. 
He  received  the  B.Eng.,  B.Eng.  (Honours)  and  M.Eng.  de¬ 
grees  (all  cum  laude)  from  the  University  of  Pretoria  in  1982, 
1983  and  1986  respectively,  and  in  1991  he  received  the  Ph.D. 
degree  from  the  University  of  Stellenbosch.  His  Ph.D.  dis¬ 
sertation  was  entitled  “Parallel  algorithms  for  electromag¬ 
netic  moment  method  formulations”.  David  is  presently  an 
Associate  Professor  of  Electrical  and  Electronic  Engineering 
at  the  University  of  Stellenbosch,  where  he  started  teaching 
in  1988  following  three  years  as  a  research  engineer  at  the 
national  research  laboratories.  His  research  interests  centre 
around  computational  electromagnetics  —  theory,  code  de¬ 
velopment  and  applications  —  and  include  the  application  of 
parallel  processing  to  this  field;  the  work  reported  here  was 
closely  linked  to  his  work  on  a  parallel  version  of  NEC2.  He 
has  worked  with  both  transputer  arrays  for  moment  method 
codes  and  the  CM-2  for  FDTD  codes.  He  recently  spent  six 
months  as  a  visiting  scholar  at  the  University  of  Arizona, 
Tucson,  AZ,  working  on  applying  computational  electromag¬ 
netics  methods  to  optical  device  modelling.  In  his  spare  time 
he  climbs  mountains,  builds  models,  board-sails,  plays  tennis 
or  squash  or  listens  to  music.  David  is  unmarried. 


39 


A  GUIDE  TO  IMPLEMENTATIONAL  ASPECTS  OF  THE  SPATIAL-DOMAIN  INTEGRAL 
EQUATION  ANALYSIS  OF  MICROSTRIP  ANTENNAS 


Louis  T.  Hildebrand  &  Derek  A.  McNamara* 
Electromagnetism  Group,  Department  of  Electrical  &  Electronic  Engineering 
University  of  Pretoria,  Pretoria,  South  Africa  0002 


Abstract:  In  the  analysis  of  microstrip  radiating 
structures  integral  equation  methods  rate  among  the 
more  accurate  approaches.  By  far  the  most  demanding 
part  in  the  use  of  these  methods  is  the  actual 
numerical  implementation.  This  paper  presents  a 
detailed  illustrated  description  of  the  numerical 
implementation  of  a  spatid  domain  mixed-potential 
integral  equation  analysis  of  microstrip  radiating 
structures. 

1.  INTRODUCTION 

In  the  analysis  of  microstrip  radiating  structures 
integral  equation  methods  rate  among  the  more 
accurate  approaches.  The  most  widely  used  is  the 
"full-wave"  electric  field  integral  equation  (EFIE). 
When  both  vector  and  scalar  potentials  are  used  in  the 
formulation  of  the  EFIE,  it  is  often  referred  to  as  the 
mixed-potential  integral  equation  (MPIE).  There  have 
been  two  paths  followed  in  the  implementation  of  the 
integral  equation  methods  to  microstrip  geometries; 
these  can  be  classified  according  to  the  manner  in 
which  the  somewhat  cumbersome  Green  functions 
appropriate  to  the  problem  are  dealt  with  numerically 
in  the  moment  method  solution  of  the  integral 
equation.  In  both  cases  one  starts  with  a  Green 
function  derived  in  analytical  form  in  the  spectral 
domain.  Expressions  for  the  moment  method  matrix 
elements  are  then  established,  consisting  of  three 
double  integrals,  over  the  spectral  domain,  domain  of 
the  testing  functions,  and  domain  of  the  expansion 
functions  for  the  unknown  current-density, 
respectively.  It  is  in  the  evaluation  of  this  expression 
that  two  different  approaches  are  used.  In  the  fiill- 
wave  spatial  domain  approach  this  is  done  by  carrying 
out  the  integration  with  respect  to  the  spectral 
variables  numerically  in  order  to  convert  the  Green 
function  to  spatial  domain  form,  and  then  proceeding 
further  with  the  other  two  integrations.  In  the  full- 
wave  spectral  domain  approach,  the  expression  for  the 
moment  method  matrix  elements  is  manipulated  in 


such  a  way  that  it  becomes  one  consisting  of  two- 
dimension^  spatial  Fourier  transforms  of  the 
expansion  and  testing  functions  (and  this  can  usually 
be  done  in  closed  form)  whose  product  with  the 
original  spectral  domain  Green  function  must  be 
integrated  in  the  spectral  domain.  In  both  cases  the 
moment  method  solution  provides  the  expansion 
coefficients  of  the  unknown  current-density  directly  in 
the  spatial  domain  of  course.  An  exposition  of  the 
difference  between  the  two  approaches  is  given  by 
Pozar  [l,Sect.IV,  Part  B]. 

Published  information  on  the  application  of  the 
integral  equation  approaches  (no  doubt  in  order  to 
satisfy  the  space  restrictions  associated  with  journal 
articles)  essentially  takes  the  form  of  some  usefiil 
suggestions  on  how  to  overcome  numerical 
difticulties,  but  has  to  stop  short  of  the  details  needed 
for  a  direct  implementation  of  the  analysis  in  a  code. 
This  is  unfortunate  since  by  far  the  most  demanding 
part  in  the  use  of  integral  equation  analyses  for 
microstrip  antenna  problems  is  the  actual  numerical 
implementation.  One  does  not  always  want  to  ask 
others  directly  for  their  computer  codes  since  these 
are  often  not  available  for  distribution  for  proprietary 
reasons,  or  the  codes  are  directly  related  to  the 
livelihood  of  the  developer  who  has  for  this  reason 
invested  a  considerable  amount  of  time  in  its 
development  (and  this  is  quite  understandable).  This 
paper  therefore  presents  a  detailed  illustrated 
description  of  the  numerical  implementation  of  the 
MPIE  analysis  of  Mosig  and  Gardiol  [2],  hereafter 
simply  referred  to  as  the  MPIE  formulation',  it  is  our 
experience  that  the  discussion  on  implementation 
given  in  [2,3]  is  the  most  detailed  available  in  the 
literature  at  present.  We  include  the  information  we 
would  have  liked  to  have  had  available  when  we 
began  implementing  the  techniques  of  [2,3].  It  will  be 
assumed  that  the  reader  has  references  [2],  [3]  and  [4] 
at  hand,  although  matters  will  be  summarised  in 
Section  2. 


‘  D.A.  McNamara  was  with  the  University  of  Pretoria.  He  is  now  with  COM  DEV  Ltd,  155  Sheldon  Drive,  Cambridge,  Ontario, 
NIR  7H6,  Canada 
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2.  BRIEF  OVERVIEW  OF  THE  INTEGRAL 
EQUATION  MODELLING  SCHEME 

Since  the  Green  functions  forming  the  kernel 
of  the  integral  equation  are  those  for  a  horizontal 
electric  current  element  (at  z=0)  in  the  presence  of  a 
grounded  dielectric  slab,  the  only  unknowns  are  the 
current/charge-densities  on  the  etched  conductors.  The 
integral  equation  is  obtained  by  "replacing"  the  etched 
conductors  by  an  equivalent  surface  current-density 
J,(x,y,0)  and  enforcing  the  boundary  condition  that  the 
total  tangential  electric  field  at  all  points  on  the  slab 
surface  occupied  by  these  conductors  be  set  equal  to 
Z,(zxJJ,  where  Z,  is  the  surface  impedance  of  the 
conducting  material.  The  validity  of  using  the  surface 
impedance  concept  has  been  carefully  studied  by 
Theron  and  Cloete  [5].  These  integral  equations  are 
expressed  in  terms  of  vector  and  scalar  potentials 
A(x,y,0)  and  $(x,y,0),  respectively,  and  thus  the 
Green  functions  us^  are  those  associated  with  these 
potentials,  namely  G^^  and  Gy,  as  described  in  [2].  A 
method  of  moments  procedure  with  sub-sectional 
expansion  and  testing  functions  is  used  to  solve  the 
MPIE,  with  overlapping  rooftop  functions  in  both  the 
X  and  y  directions  selected  for  the  expansion  functions 
of  the  surface  current-density.  These  are  placed  as 
follows:  the  surface  of  the  etched  conductor  is  divided 
into  elementary  charge-density  cells,  with  two  adjacent 
charge-density  cells  forming  a  current-density  cell. 
Figure  1(a)  shows  two  charge-density  cells.  The 
current-density  is  represented  by  overlapping  rooftop 
expansion  functions.  So-called  razor  test  functions, 
shown  in  Figure  1(b),  which  extend  along  segments 
linking  the  centres  of  adjacent  charge-density  cells, 
are  used.  Although  not  a  fondamental  limitation  of  the 
technique,  we  have  here  arranged  matters  so  that  the 
antenna  structure  must  be  segmented  into  rectangular 
cells  of  equal  size;  this  simply  eases  some  of  the 
already  substantial  computation^  burden. 

3.  THE  DISCRETISED  (MATRIX)  FORM  OF 
THE  INTEGRAL  EQUATION 

Let  there  be  M  of  the  x-directed  rooftop  functions, 
and  N  of  the  y-directed  ones.  Their  introduction, 
along  with  the  razor  testing  functions  extending  from 
the  centres  of  two  adjacent  charge  cells  (that  is,  along 
C,^  in  Figure  1),  converts  the  integral  equation  to  the 
matrix  equation  [C][I]  =  [V]  which  is  written  [3,4]  in 
terms  of  sub-matrices  as 


c  C’^ 

h 

1 

Vw' 

’  J^o 

vy 

(a)  Basis  fiinctioii 


(b)  Razor  testing  fanction 


Figure  1  X-directed  current  cell  with  a  rooftop  basis ftmction  and 
razor  testing  ftmction.  The  associated  charge  distribution  over  the 
current  cell  is  also  shown. 


Here  the  C”,  C^,  and  C”'  are  sub-matrices  of  size 
MxM,  MxN,  NxM  and  NxN,  respectively.  Thus 
matrix  [C]  is  of  size  (M-l-N)x(M-l-N).  The  elements 
of  the  sub-matrices  are  given  by  [3,4] 

r^(r*lr;^)  -  r^(r>;) 

+  r^(r,>;,)  +  rp(r,>;y)  ] 

z  ® 

i  -  I...M  J  -  1...M 


-  —  — -  [  -  r,,(r,,lrj,y)  -  r,,(r,jlryy) 

+  r,,(r;,lr;p  +  TyirJrjj)] 
i  -  \  ...M  ,}  -  I...N 


(3) 


In  the  above,  Zp  is  the  free  space  characteristic 
impedance,  kp  the  free  space  wavenumber  and  6^  the 
Kronecker  delta.  The  Z,  is  a  surface  impedance 
accounting  for  the  finite  conductivity  as  well  as 
surface  roughness  and  finite  thickness  of  the  metallic 
upper  conductor  [3].  An  expression  for  is  obtained 
by  interchanging  the  couplets  (x,y),  (a,b)  and  (M,N) 
within  the  expression  for  C^;  reciprocity  requires  that 
for  cells  of  equal  size  C]]  =  C^. 

P  and  P  are  column  sub-vectors  of  size  Mxl 
and  Nxl,  respectively.  Their  elements  are  the 
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coefficients  of  the  x-directed  and  y-directed  expansion 
functions  for  the  current-  density.  The  excitation 
column  sub-vectors  and  are  similarly  of  size 
Mxl  and  Nxl,  respectively.  Thus  the  column  vectors 
[I]  and  [V^®^]  are  both  of  size  (M+N)xl. 

The  precise  expression  for  the  elements  of  [V] 
depends  on  the  type  of  excitation  geometry  being 
considered.  Here  a  coaxial  probe  feed  will  be 
considered  for  which  a  simple,  yet  sufficientiy 
accurate  model  was  introduced  in  [3].  This  model 
assumes  that  the  current  on  the  coaxial  probe  is 
constant  and  is  therefore  only  accurate  for  thin 
substrates  (up  to  about  X/10  thick).  According  to  this 
model  the  excitation  current  spreads  over  a  single 
charge  cell  as  described  by  the  following  expression: 

\  “  I  ,4) 

More  complex  feed  models  valid  for  thick  substrates 
are  discussed  in  [6].  The  associated  excitation  surface 
charge-density  distribution  over  the  charge  cell  is 
given  by  a  rectangular  pulse  of  value  I/jwab  where  I 
is  the  total  current  carried  by  the  inner  coaxial 
conductor.  Figure  2  gives  an  illustration  of  the  electric 
surface  current-density  and  charge-density 
distributions  associated  with  the  coaxial  probe  feed 
model;  note  that  only  the  x-directed  component  is 
shown  in  the  figure  for  the  sake  of  clarity.  It  has  been 
found  [7]  that  the  contribution  of  the  excitation 
current  (as  opposed  to  the  charge)  to  can  usually 
be  neglected;  otherwise  it  must  be  computed  from  the 
expression  for  T”  (T^^O  given  in  (8)  with  (J^) 
replacing  T,  (Ty).  If  we  recognise  that  the  matrix 
elements  according  to  (3)  represent  the  effect  of  a 
charge  doublet  (Figure  1)  integrated  along  a  test 
segment,  then  since  with  the  above  coaxial  probe 
model  the  excitation  is  a  single  pulse  of  charge 
(Figure  2),  the  elements  of  the  excitation  vector  may 
be  seen  to  be  approximated  by  [7] 

(5) 


where  feedpoint  r  is  located  at  the  centre  of  a  charge 
cell.  The  moment  method  matrix  C  is  ill-conditioned 
due  to  the  fact  that  some  rows  are  almost  linear 
combinations  of  three  other  rows.  Therefore  careful 
evaluation  of  its  elements  is  necessary  to  obtain  results 
of  sufficient  accuracy.  On  the  other  hand,  C  is 
diagonally  dominant,  Aerefore  less  stringent  accuracy 
requirements  apply  to  the  off-diagonal  elements. 
Certain  approximations  may  then  be  considered.  One 
of  these  approximations  involves  the  integral  term  in 
(2)  which  may  be  approximated  as  follows; 

(6) 

The  validity  of  this  approximation  is  illustrated  in 
Figures  3  and  4  showing  the  real  and  imaginary  parts, 
respectively,  of  both  terms  in  (6).  In  fact,  this 
approximation  may  even  be  considered  for  the 
diagonal  elements,  since  the  contribution  of  the  vector 
potential  to  the  value  of  the  matrix  element  is 
overshadowed  by  that  of  the  scalar  potential  [7].  A 
conjugate  gradient  method  [8]  is  used  for  the  iterative 
solution  of  the  moment  method  matrix  equation. 

In  (2),  (3)  and  (5)  F”  and  Fy  are  discrete  Green 
functions  (not  to  be  confiised  with  the  Green  functions 
themselves)  which  have  as  sources  complete  basis 
functions  as  opposed  to  conventional  elementary  point 
sources.  Fi’‘(r|rj)  is  then  the  x-component  of  the 
vector  potential  at  r  created  by  an  x-directed  rooftop 
distribution  of  surface  current  at  Tj,  whereas  Fv(r  |  r,j) 
is  the  scalar  potential  at  the  same  observation  point 
resulting  from  a  rectangular  distribution  of  unit 
surface  charge  at  r,j.  The  discrete  Green  functions  are 
now  defined  by  the  following  dimensionless 
expressions  [3]: 

r^(rlr^)-/  ^G^(rlr')n(r'-r^)*^dS'  (7) 

where  and  €„  are  the  permeability  and  permittivity 
of  free  space.  A  similar  expression  to  that  given  in  (8) 
holds  for  Fi’'.  S,^  represents  the  surface  of  the  current 
source  cell  centred  at  (Figure  1)  while  the  charge 
source  cell  centred  at  extends  over  S<^.  11  is  a  two- 
dimensional  unit  pulse  function  over  and 
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Tj(r)=l  — |x|/a  over  S^.  The  Green  functions 
appearing  in  the  above  expressions  for  the  discrete 
Green  fiinctions,  are  given  by  the  following 
expressions: 

m 

G„(rlr') - —  f  UXR)  ^  ^  W 

^  27r€oi  ° 


GT(r\r')  -  G»(rlr')  -  ^ / J^(kR)-^e-'^^dX 

(10) 

where  Jq  is  the  Bessel  function  of  zero’th  order  and 
first  kind,  DTE=Mo+wcoth(M/i),  DjM=e^+tA3i\h(uh) 
and  N=Mo+Mtanh(Mh);  Uo=(X^-k^y'^  u=(y-e}i^y'^, 
h  is  the  substrate  thickness  and  R=  |r-r'|.  For  the 
case  of  a  lossy  dielectric,  er=e/(l  — jtan6)  where  tanS 
is  the  loss  tangent.  Figure  5  illustrates  the  effect  of  the 
distance  between  source  and  observer  (R)  on  the 
Green  fiinction  integrands;  the  real  part  of  the 
normalised  G“  integrand  is  shown  for  R=0.5Xo  and 
O.OSXq.  It  can  be  seen  that  when  the  observer  (r) 
approaches  the  source  point  (r'),  the  zero’s  of  the 
integrand  move  away  from  the  origin.  In  the  limiting 
case  (R=  |r-r'|^)  the  first  zero  tends  to  infinity, 
producing  a  non-oscillating,  non-zero  function  to  be 
integrated  over  a  semi-infinite  interval.  This  produces 
the  singularities  in  the  Green  functions  at  the  source 
point.  Now,  in  the  evaluation  of  the  discrete  Green 
functions  the  situation  arises  where  the  observation 
point  r  for  the  potential  due  to  a  particular  source  cell 
falls  within  its  own  source  cell  boundaries.  For  this 
reason  the  integrands  for  the  discrete  Green  functions 
are  singular  at  r  due  to  the  Green  functions  becoming 
singular.  This  situation  is  illustrated  in  Figures  6  and 
7.  We  will  call  any  evaluation  of  the  potential  due  to 
a  given  source  cell  at  any  point  in  its  own  source  cell 
a  "self  term".  By  now  writing  (7),  in  the  evaluation  of 
the  self  term,  as 
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it  is  possible  to  evaluate  the  scalar  potential  discrete 
Green  function  self  term  since  the  difference  term  in 
(11)  is  numerically  integrable  as  illustrated  in 
Figures  8  and  9.  Even  though  the  difference  term  is 
well-behaved  at  r'=r,j  both  terms  comprising  the 
difference  term,  i.e.  Gy  and  Gy.,  become  singular. 
The  evaluation  of  the  difference  term  at  the  source 
point  is  therefore  done  as  follows: 
Gv(rojk,4+fi)  -  Gv.(r^k.j+6)  where  6-*0.  In  this 
way  an  accurate  estimate  of  the  well-behaved 
difference  term  may  be  obtained  at  the  source  point. 
In  similar  fashion  the  singularity  in  the  real  part  (the 
imaginary  part  of  the  F”  integrand  does  not  exhibit 
singular  behaviour  even  in  self  term  evaluations)  of 
the  vector  potential  discrete  Green  function  may  be 
dealt  with.  In  this  case,  however,  an  expression  for 
FjJJ  could  not  be  found  in  the  literature,  however, 
analytical  integration  yields: 


Tt  2n  ^2  4) 


“(■T 


2  2 

-^^Intanf— + — )+ -^^^fcosecfa  J -cosec(a  i)! 
4na  2  4  16na^  •* 
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a, 
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2na  2  4 

(13) 


where  ai=tan“^[b/(2(a-r,))],  a2=ta°~^[b/(2(a+rj)] 
and  a3=tan“‘[b/(2r^].  IThis  expression  was  derived 
for  an  x-directed  current  cell  centred  at  (0,0)  with  an 
observer  on  the  source  cell  at  Therefore, 
numerical  techniques  required  in  the  self  term 
evaluation  of  both  the  discrete  Green  functions  have 
now  been  discussed.  For  the  off-diagonal  matrix 
elements  approximations  for  the  discrete  Green 
functions  were  suggested  by  Mosig  and  Gardiol  [3]. 
It  is  worth  mentioning  that  Aese  approximations  have 
also  been  implemented  successfully  by  the  present 
authors. 


Figure  2  Electric  surface  current  and  charge  distributions 
associated  with  the  coaxial  probe  feed  model. 
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Figure  3  Real  parts  of  the  actual  integral  ( - )  and 

cq/projdmation  (6)  ( — O — 0 — )  suggested  by  Mosig  and 
Gardiol  [1];  f=1.206  GHz,  e=4.34-j0.0868;  h=0.8mm  and 
a=b— 6.666  mm. 
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1  3  5  7  9  x/» 


Figure  4  imaginary  parts  of  the  actual  integred  ( - ) 

and  approximation  (6)  ( — O — 0 — )  with  the  corresponding 
real  part  shown  in  Figure  3. 
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Figure  5  Normalised  integrand  of  Of‘(r/r')forR/hf,=0.5( - ) 

and  R/Ko=0.05(—X—y.—)  where  h/i^-0. 07  and f- 1.206  GHz 


Figure  6  Real  part  of  the  integrand,  brfore  pole  extraction, 
in  the  expression  for  the  Fy  selfterms. 


Figure  7  imaginary  part  of  the  Fy  integrand  corresponding 
to  Figure  6.  Note  that  the  maxitnum  value  is  infinite. 


Figure  8  Real  part  of  the  difference  term  integrand  in  the 
evaluation  of  the  Fy  selfierm. 


Figure  9  imaginary  part  of  the  difference  term  integrand 
corresponding  to  the  real  part  shown  in  Figure  8. 

4.  DETAILED  ILLUSTRATION  AND 

DISCUSSION  OF  THE  NUMERICAL 
EVALUATION  OF  THE  ASSOCIATED 
INTEGRALS 

We  next  consider  the  evaluation  of  the  Green 
functions.  Figure  10  shows  the  integrand  of 
(2Trlfio)GX^  as  a  function  of  X/ko  at  1.206  GHz  for 
e,=4.34— j0.0868,  and  the  particular  circumstances 
h/Xo=0.07  and  R/Xo=0.5  (Xq  is  the  free  space 
wavelength).  The  scalar  potential  counterpart,  i.e.  the 
integrand  of  2‘reoGv  as  a  function  of  X/ko,  is  shown  in 


Figure  11.  As  seen  from  these  figures,  discontinuous 
derivatives,  singularities  and  oscillations  in  the  Green 
function  integrands  pose  distinct  problems, 
complicating  the  use  of  numerical  integration 
techniques.  Mosig  and  Gardiol  [9]  therefore  suggested 
that  the  semi-infinite  integration  intervals  in  (9)  and 
(10)  be  subdivided,  allowing  these  problems  to  be 
addressed  separately.  For  convenience,  the  Green 
function  integrals  are  written  as 

*0 

-  f  F(X)dX  + 

0  0 

os 

f  F(X)dX+  f  F(X)dX 

(14) 

The  numerical  difficulties  in  the  evaluation  of  the 
Green  function  integrands  in  each  of  these  sub¬ 
intervals  will  now  be  discussed,  and  the  use  of 
proposed  solutions  illustrated. 

4.1  Interval  XE[0,ko] 

The  term  «o=(X^-k^''^  appearing  in  the 
expressions  for  Dte,  and  N  introduces  branch 
points  in  the  Green  function  integrands  at  X=ko. 
Manifestations  thereof  are  the  discontinuities  in  the 
derivatives  seen  at  X=ko  (Figures  10  and  11). 
Standard  numerical  integration  routines  may  be 
inefficient  in  the  integration  of  functions  with  such 
discontinuous  derivatives.  To  obtain  accurate 
numerically  integrated  estimates  for  the  Green 
functions  over  the  first  interval,  Mosig  and  Gardiol 
[4]  have  proposed  the  substitution  X=kocos  t.  Suppose 
the  integrands  may  be  written  as  F(X),  then  this 
substitution  implies  that 

*6  2k 

J^F(X)dX‘‘  J  F{kQCOst)(-kQsmt)dt 

0  3kI2 

(15) 

Figures  12  and  13  show  F(koCOS  t)(-koSin  t)  as 
fimctions  of  t,  with  F  representing  the  integrands  of 
(27r/;io)GA*  and  2ireoGv  shown  in  Figures  10  and  11 
respectively.  The  integrands  are  found  to  be  smooth 


f  F(X)dX 
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and  easily  integrable.  In  this  way  then,  it  is  possible 
to  deal  with  the  effect  of  the  branch  points  in  the 
Green  function  integrands  at  X  =  ko. 

4.2  Interval  X€[ko,ko(e/)*^ 

Since  the  integrand  of  Gy  (9)  contains  Djm  ui  the 
denominator,  a  singularity  due  to  the  existence  of  the 
dominant  TMq  surface  wave  mode,  j^pears  in  the 
interval  [ko  ,  ko(e,')‘'^],  as  seen  from  Figure  11. 
According  to  a  pole  extraction  technique  described  in 
[4],  the  2x-6oGv  integrand  may  be  expanded  as 
follows: 

UXR)XN 

. - 

"re  "tv 

where  F,i,g(X)  =  Res/[X  -  (Xp  +  j  i^p)] ;  Res  is  the  residue 
of  F(kp)  (where  kp=X  +  iv)  at  pole  kpp=Xp  +  ji^p 
(F(X)  is  a  special  case  of  FO^J).  This  leaves  k^p  and 
Res  to  be  determined.  The  present  authors  have  used 
Muller’s  method  [10]  with  deflation  to  determine  k^p  - 
a  root  of  the  univariate  complex  function  DTM(kp). 
This  method  is  numerically  implemented  in  the  IMSL 
routine  ZANLY  [11]. 

The  residue  may  furthermore  be  determined  as 
follows  : 

(16) 

where  kp(t)=X(t)+j»'(t),  t  G  [a,b]  and  C  is  any  closed 
path  around  the  pole  at  k^p.  A  condition,  however,  is 
that  function  F(kp)  must  be  analytic  inside  C  except  at 
kpp  [12].  Since  C  may  be  represented  by  any  closed 
path  around  k^p,  it  is  convenient  to  choose  a  circle. 
Then  we  have  X(t)=rcos(t)  and  »'(t)=rsin(t)  as  the 
parametric  equations;  t  G  [0,27r)  and  r  the  radius  of 
the  circle.  The  value  of  r  is  not  important,  provided 
F(kp)  is  analytic  inside  the  borders  of  the  circle  except 
at  kpp.  Figure  14  shows  the  normalised  distance 
I  kpp/1^  —  1 1  between  the  pole  at  k^  and  the  branch 
point  at  X=ko  as  a  function  of  dielectric  thickness 
(h/Xo).  It  can  be  seen  that  in  the  case  of  electrically 
thin  substrates,  the  pole  due  to  the  surface  wave  is 
very  close  to  the  branch  point.  In  such  cases,  radius 
r  must  be  chosen  carefully  to  avoid  the  inclusion  of 
the  branch  point  into  the  borders  of  the  circle 
inclusion  thereof  violating  the  analyticity  of  F(kp). 
With  the  use  of  (16)  and  careful  selection  of  r,  it  is 


now  possible  to  numerically  determine  a  value  for 
Res.  For  instance,  at  1.206  GHz  (ko=25.2753)  for 
h/Xo=0.07,  R/Xo=0.5,  e,=4.34-j0.0868  and 

r=  1.0153  m-‘  we  have  Res = 0.47323 -jO.Ol 8 15 
where  kpp=27.3059-j0.052039.  Once  k^p  and  the 
residue  are  known,  the  pole  in  the  Gy  integrand  may 
be  extracted  according  to  the  technique  describe 
above;  that  is 


f  F(X)dX-  f  IFW-F^mV)^ 

*b  *b 

*bv^ 


F,i„g  is  analytically  integrable  [4]  while  the  difference 
term,  which  is  a  well-behav^  function,  may  be 
integrated  numerically.  In  the  lossless  case  (i.e.  for 
tan6=0)  we  have  kpp=Xp,  in  other  words,  the  pole  due 
to  Dtm  lies  on  Ae  X-axis  which  is  the  path  of 
integration.  Therefore,  at  Xp  both  terms  constituting 
the  difference  term  in  (17),  become  singular. 
However,  by  evaluating  F(X+6)  -  F.i„g(X-t-6)  where 
5  -►  0  it  is  possible  to  obtain  an  accurate  estimate  of 
the  well-behaved  difference  term  at  Xp.  In  the  case  of 
lossy  dielectrics,  the  pole  at  k^p  does  not  lie  on  the 
path  of  integration  (Figure  1 1),  nevertheless,  strong 
variations  in  the  integrand  due  to  the  pole  still  require 
application  of  this  pole  extraction  technique.  An 
infinite  derivative  in  the  difference  term  integrand  at 
X=k<,  may  then  be  eliminated  by  the  substitution 
X=kocosh  t.  The  real  and  imaginary  parts  of  the  Gy 
difference  term  integrand,  after  substitution,  are 
shown  in  Figure  15;  the  singularity  which  has  been 
extracted  is  visible  in  Figure  1 1 .  If  the  integrand  does 
not  contain  in  the  denominator,  the  substitution 
may  nevertheless  be  performed  to  obtain  a  smooth 
integrand  at  X=ko.  This  is  the  case  for  the  Green 
fiinction  Gi*(r|r'),  for  which  Figure  16  gives  an 
illustration  of  the  integrand  shown  in  Figure  10,  after 
substitution.  At  this  point,  we  are  able  to  accurately 
determine  the  integrals  in  (9)  and  (10)  for 
0  ^X^  ko(e/)''^  This  leaves  the  interval  X>ki,(c/)''^ 
which  is  the  subject  of  the  following  section. 
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4.3  The  method  of  averages  in  the  interval 
xe[ko(6/)^,oo] 

As  seen  from  Figures  10  and  11  the  Green  function 
integrands  show  oscillatory  behaviour  in  this  interval. 
Furthermore,  these  integrands  have  envelopes  which 
converge  very  slowly  and  therefore  standard 
integration  routines  (such  as  the  trapezium  rule  and 
Gauss  quadrature)  prove  to  be  very  inefficient  since  a 
large  number  of  integration  points  is  required  to 
achieve  reasonable  accuracy.  Mosig  and  Gardiol  [9] 
have  found  a  technique  known  as  the  method  of 
averages,  introduced  by  Hurwitz  and  Zweifel  [13],  to 
be  suitable  for  application  to  Sommerfeld  integrals 
appearing  in  microstrip  problems.  This  method  is 
based  on  the  decomposition 

fg(.KR)fa)d^  -  Y,  f  gaR)fa)di 

a  "“O  a*Hpli 

(18) 

where  g(|R)  is  an  oscillating  function  with  period  p 
and  f(^)  a  smooth,  non-oscillating  function  which 
behaves  asymptotically  as  0(|“);  the  integrand 
therefore  diverges  for  a>0.  Although  Bessel 
functions  of  the  first  kind  -  appearing  in  the  Green 
function  expressions  -  are  not  strictly  periodic,  the 
method  of  averages  may  still  be  applied.  However, 
since  the  zero’s  of  Bessel  functions  are  not  known  off¬ 
hand  the  large-argument  approximation 
J„C^)»[2/(7rXR)]''^cos(XR-ir/4-n7r/2)  is  used  to 
estimate  their  zero’s.  Then  we  have 


(«-l)  +  0.75  +  - 
2 


11 

R 


(19) 


where  approximates  the  m’th  zero  of  J„(XR). 
Table  3.1  gives  an  indication  of  the  accuracy  of  this 
approximation  in  comparing  (m=  1, 2,3,4)  with  the 
actual  zero’s  of  Jo(X)  determined  numerically  with  the 
IMSL  routine  ZREAL  [11].  A  question  about  the 
validity  of  this  approximation  for  small  R  values 
(source  and  observation  points  close  to  each  other) 
might  well  arise  at  this  point.  In  order  to  answer  this 
question,  definite  integrals  for  which  the  answers  are 
known,  were  evaluated  by  the  authors  using  (19)  in 
the  method  of  averages.  We  know  from  [14]  that 


{  ’  « 

for  n  >  -1.  This  identity  was  confirmed  for  R  values 
ranging  from  2.0  x  10"*  to  2.0  and  for  both  n=0  and 
n=l.  Further  tests  performed  by  the  authors  on 
similar  identities  lead  to  the  conclusion  that  (19)  will 
not  introduce  a  significant  error,  even  for  small  values 
of  R,  when  used  in  the  method  of  averages. 


m 

m’th  zero  of 
Jo(X) 

km 

error 

(%) 

1 

2.4(H825 

2.356194 

2.02 

2 

5.520078 

5.497787 

0.40 

3 

8.653727 

8.639379 

0.16 

4 

11.79153 

11.78097 

0.08 

Table  3.1  Zero’s  of  Jo(X)  determined  with  ZREAL 
[11]  and  (19)  respectively. 


Figure  17  gives  the  method  of  averages  in  the  form  of 
a  flowchart.  The  first  step  is  to  perform  an  integration 
over  a  half  cycle  to  determine  1^. 


(m-l) 

a 

iLi*  /  sam^di  ('«>i) 


(21) 

where  is  the  m’th  successive  zero  of  the  oscillating 
function  g,  with  >  a.  I^_i  is  then  determined 
through  the  use  of  a  weighted  mean 


1  rl 


+ 


w. 


w-l 
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(22) 
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with  bofli  li  and  Ii_,  having  been  determined 
previously  through  integration  and  with  the  weights 
given  by  wj^  =  (li/|  In  general,  (22)  is  given 

by 


t-l  r*-l 

*  W_  i_  + 


w*'*  /*•* 
^■1+1  ^«+l 


w‘-‘  . 


w. 


.n 


(23) 


■141 


If  m-1  1  then  I’_2  is  calculated  through  the  use  of 

(23)  and  the  process  is  repeated  until  we  have  Ij 
(where  k=m)  which  will  give  an  approximation  of  the 
actual  value,  I(R),  despite  the  fact  that  no  additional 
integrations  were  performed.  If  the  error  criterion  is 
not  met,  li  (where  m  =  k+1)  may  be  determined 
through  integration,  (21),  and  Ij^'  through  repeated 
implication  of  (23).  In  this  way,  an  estimate  for  I(R) 
may  be  obtained. 

In  the  determination  of  the  discrete  Green  function 
self  terms,  the  situation  arises  where  the  distance 
between  source  and  observer  (R)  tends  to  zero.  The 
effect  of  this,  as  discussed  earlier  and  illustrated  by 
(19),  is  that  the  zero’s  of  the  oscillating  Green 
function  integrand  move  further  out  from  the  origin 
along  the  X-axis.  In  the  method  of  averages,  weights 
are  determined  according  to  the  zero’s  of  the 
oscillating  function,  i.e.  Now  it  is 

imparent  that,  for  ^„  >  >  and  k  a  positive  integer 
which  may  be  large,  the  weight  w‘  could  become  a 
very  large  real  number,  creating  possible  numerical 
difficulties  (numeric  overflow  or  round-off  errors,  for 
instance).  The  present  authors  have  addressed  this 
problem  through  mathematical  manipulation  of  (23) 
after  which 


jk-l 
jk  •‘m 
■*in  ” 


up 


up 


,k-l 
*in  +  l 

rk-l 


■JB 


/J 


(24) 


where  /3=  10exp{(o(  +  2-k)[logio(li/{„+i)- 
logio(^i/^  J]}  is  now  a  manageable  real  number  even 
though  w^“'  and  w^^:;:}  are  very  large.  Consider  for 
example  a  situation  where  ^,=200,  ^30=8x10^, 
^31=8.1  xK)*,  k=50  and  a=0;  then  we  have 
wg=7.923xl0’«  and  wj?=  1.438x10^  whilst  /S  = 
1.815  !.  Therefore,  to  avoid  numerical  difficulties 
which  may  be  encountered  in  a  straightforward 
application  of  (23),  we  propose  (24)  as  a  means  to 
determine  Ij^. 


Figure  10  Real  ( - )  and  imaginary  ( — a — a — )  integrands 

of  (Ix/n^G^for  e=4.34-j0.086S,  h/ko=0.07  and  R/\o=0.5. 


GHz. 


Figure  12  Real  ( - )  and  imaginary  (• — a — a — )  integrands 

cf  (2t/iio)G"  for  X  €  [0,kj  after  substitution  \=k(fios  t. 
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Figure  13  Real  ( - )  and  imaginary  (—a — a — )  integrands 

oflTCgGy  in  the  interval  X  €  [0,kj  after  substitution  'k=kgCos  t. 


Figure  16  integrand  of  (2ie/no)Gf  (real  (- 
imaginary  ( — a — a — a — ))  after  substitution. 


-)  and 


Figure  14  Normalised  distance  between  the  pole  due  to  Djy  and 
the  branch  point  due  to  Ug  as  a  junction  of  dielectric  thickness,  i.e. 

I  I  (h/kc). 


imaginary  ( — a — a — a — ))  after  substitution  'K=kgcosh  t. 

5.  INTERPOLATION 

As  discussed  in  the  previous  section,  the  evaluation 
of  the  Green  functions  require  implementation  of 
numerical  integration  routines.  This  is  in  addition  to 


Figure  17  Flowchart  for  the  method  of  averages 

the  numerical  integration  which  has  to  be  performed 
in  the  evaluation  of  the  discrete  Green  functions 
(which  contain  the  Green  functions).  With  these 
discrete  Green  functions  it  is  now  possible  to  compute 
the  moment  method  matrix  elements  whose  solution 
yields  the  unknown  coefficients  !,( and  ly.  However,  in 
order  to  reduce  the  computation  time,  Mosig  and 
Gardiol,  noting  the  fact  that  the  relevant  Green 
functions  are  only  dependent  on  the  distance  between 
source  and  observer  (and  not  the  relative  orientation), 
consequently  suggest  the  use  of  an  interpolation  table 
to  ev^uate  the  Green  functions.  This  interpolation 
table  consists  of  a  discrete  set  of  Green  function 
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values  tabulated  against  a  set  of  distances  R;  where 
i=l..P  and  <  R^  <  R^.  Interpolation  is  then 
used  to  evaluate  these  functions  at  any  distance 
ranging  from  R^  to  R^  (the  maximum  linear 
dimension  of  the  structure).  It  is  suggested  that 
convergence  tests  be  performed  to  determine  optimum 
values  for  P  (typically  from  50  to  250)  related  to 
specific  problems.  Since  the  value  of  R^  has  an 
effect  on  the  self  term  evaluation  of  the  discrete  Green 
functions,  care  should  be  taken  when  choosing  values 
for  R^.  However,  numerical  experiments  performed 
by  the  present  authors,  showed  differences  on  the 
order  of  0.1%  in  the  self  term  values  computed  for 
R^=10~*  and  for  R^=10  “.  Furthermore,  careful 
consideration  should  be  given  to  the  computation  of 
Ri.  In  this  regard,  the  following  expressions  have  been 
used  to  good  effect: 


is  shown  in  Figure  18  which  illustrates  the  r^id  sign 
changes  and  reduced  integration  interval  just  spoken 
of.  From  this  figure  it  is  clear  that  standard  numerical 
integration  routines  would  not  be  capable  of  yielding 
accurate  estimates  of  the  Green  fiinction  integrals  with 
large  z  values.  We  therefore  resort  to  asymptotic 
techniques  to  obtain  approximate  analytic  solutions. 
Application  of  the  method  of  steepest  descent  [15]  to 
these  integrals  yields  the  following  far-field  radiation 
expressions  for  arbitrarily  shaped  etched  antennas  in 
terms  of  the  coefficients  I,  and  ly  [3]: 

«e-Gl*(rlO)i:fl/y‘^"-GlVlO)i:W,,e''^^  W 

i-i  j-i 

i-i  j-i 

where  gk=XkSin0cos(^+ykSin0sin<^  (k=ij);  Xk  y^ 
are  the  centre  coordinates  of  the  k’th  current  cell.  We 
also  have 


for  R^<R<Ri^  and 


(25b) 

for  Ri^<R<R,^  where  Ri„ter=(a^+b^‘'^.  With 
interpolation  it  is  possible  to  reduce  the  computation 
time  without  sacrificing  significant  accuracy. 


6.  COMPUTATION  OF  FAR-ZONE 
RADIATION 

In  solving  for  the  surface  current  distribution,  we 
are  interested  in  the  case  where  both  source  and 
observer  are  located  on  the  air-dielectric  interface  i.e. 
where  z=0.  For  far-field  radiation  computation,  on 
the  other  hand,  this  will  not  be  the  case.  If  the 
radiator  is  placed  in  the  xy-plane,  far-field  radiation 
on  broadside  now  implies  that  z  becomes  large  and 
then  the  following  situation  arises:  since  Uq  is  purely 
imaginary  in  the  interval  0  ^  X  <  ko,  and  z  a  very  large 
real  number,  exp(— j(k2-X^*'^)  causes  rapid  sign 
changes  in  this  interval.  Furthermore,  since  Uq  is  real 
for  X^kfl  we  have  exp(-UoZ>^  and  thus  the 
integration  interval  in  effect  reduces  to  0^X<ko.  An 
example  of  a  Green  function  integrand  for  z=5.0  m 
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T  -  ,/(e,-sin2e) 


Figure  18 

Figure  18  Normalised  real  integrands  of  O'/  (  i.e. 
JgftJi)\e3p(-u^/Dj^for  R/kg=0.5  and  h/Xg=0. 07 at f  =1.206  GHz 
and  e,=4.34-j0.0868  )for  z  =  0.0  and  5.0  m. 

7.  SUMMARY 

The  MPIE  formulation  rates  among  the  accurate 
integral  equation  analysis  techniques  for  microstrip 
radiating  structures.  After  a  brief  overview  of  the 
modelling  scheme,  numerical  techniques  used  by  the 
authors  in  the  implementation  of  the  formulation  were 
discussed:  a  detailed  and  illustrated  discussion  of  the 
numerical  evaluation  of  the  required  Green  functions 
was  given,  an  interpolation  method  used  to  improve 
the  computational  efficiency  was  discussed  and  finally 
the  computation  of  far-field  radiation  patterns  was 
considered.  A  FORTRAN  computer  code 
implementation  of  the  MPIE  formulation  is  available 
on  request  from  the  second  author. 
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Abstract^ 

A  characteristic-based,  windward®  numerical  pro¬ 
cedure  for  solving  three-dimensional  Maxwell  equations 
in  the  time  domain  has  been  successfully  ported  to  the 
Intel  Touchstone  Delta  multicomputer.  The  numerical 
results  by  concurrent  computation  duplicated  the  earlier 
simulations  of  an  oscillating  electric  dipole  on  a  vector 
processor  and  compared  well  with  the  exact  solutions. 
The  parallelized  code  is  scalable  up  to  512  nodes  and 
incurs  only  up  to  7.6%  performance  degradation.  The 
sustained  data  processing  rate  is  clocked  at  6.551  Gi- 
gaops.  However,  the  data  I/O  process  is  unscalable  on 
the  shared  memory  system. 

Nomenclature 


E 

Electric  field  intensity 

H 

Magnetic  field  intensity 

hj,k 

Indices  of  discretized  grids 

J 

Electric  current  vector 

L 

One-dimensional  difference  operator 

n 

Index  of  time  level 

t 

Time 

W 

One-dimensional  characteristic 

x,y,z 

Cartesian  coordinates 

r,e,4> 

Spherical  coordinates 

€ 

Electric  permittivity 

Magnetic  permeability 

A 

Eigenvalue 

1  Introduction 

The  improvement  of  numerical  efficiency  is  one 
of  the  urgent  needs  of  computational  electromagnetics 
(CEM)  in  aircraft  signature  technology.  In  this  area  of 
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^The  windward  difference  approximation  is  achieved  by  using 
the  one-sided  stencil  for  computing  the  difference  quotient  accord¬ 
ing  to  the  sign  of  the  coefficient. 


applications,  the  CEM  simulations  are  generally  more 
computationally  intensive  than  computational  fluid  dy¬ 
namics  (CFD)  problems.  A  part  of  the  reason  is  that 
the  numerical  accuracy  requirement  for  CEM  tends  to 
be  more  stringent  than  for  CFD.  For  instance,  a  desir¬ 
able  predictive  dynamic  range  can  be  as  high  as  60  dB 
over  broad  viewing  ranges  [1] .  Typically  for  wave  prop¬ 
agation,  a  suitable  numerical  resolution  requires  each 
wavelength  to  be  supported  by  at  least  an  order  of  ten 
grid  points  or  more.  Thus,  for  a  scatterer  with  high 
refraction  indices  and  material  complexities,  this  more 
than  ten  grid  points  per  wavelength  requirement  trans¬ 
lates  into  the  need  of  a  computing  system  with  astro- 
nomic.al  computing  speed  and  memory  size. 

The  heavy  demand  on  computer  speed  and  mem¬ 
ory  for  CEM  simulations  can  be  illustrated  by  the  fol¬ 
lowing  example.  For  a  fixed  incident  angle,  the  numer¬ 
ically  generated  signature  of  a  modern  fighter  config¬ 
uration  at  1  GHz  requires  approximately  fifty  million 
grid  points®  to  produce  a  ten  grid  points  per  wavelength 
resolution.  For  each  grid  point,  the  values  of  at  least 
three  coordinates,  six  field  components,  and  nine  direc¬ 
tion  cosines  of  a  general  curvilinear  coordinate  system 
need  to  be  stored  [2, 3, 4,5].  This  results  in  a  total  of 
nearly  one  billion  memory  allocations.  A  typical  CEM 
code  operating  on  a  100  Megaops  (10®  floating  point 
and  integer  operations)  single  vector  processor  can  pro¬ 
cess  data  at  approximately  a  rate  of  three-tenths  of  a 
microsecond  per  grid  point  per  time  step.  At  this  rate, 
a  fighter  configuration  will  require  almost  five  hours  of 
computing  just  to  advance  the  solution  to  a  new  time 
level.  In  order  to  complete  a  signature  simulation,  usu¬ 
ally  multiple  look  angles  and  hundreds  of  time  steps 
are  required.  As  a  result,  the  large  computer  memory 
requirement  and  solution  time  has  rendered  the  use  of 
conventional  data  processors  for  solving  CEM  problems 
impractical. 

The  goal  of  this  paper  is  to  provide  an  efficient 
way  of  solving  the  3D  time  domain  Maxwell’s  equations 
in  differential  form  by  combining  novel  numerical  algo- 
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rithms  and  concurrent  computing  technology.  In  the 
recent  years,  the  CFD  community  has  made  significant 
progress  in  the  area  of  algorithm  development  [6,7,8]. 
Numerical  algorithms  for  solving  hyperbolic  equations^ 
from  the  CFD  discipline  have  been  adopted  for  solving 
three-dimensional  Maxwell  equations  in  the  time  do¬ 
main  [2, 3, 4, 5].  Among  these,  the  characteristic-based 
algorithm  is  found  to  be  most  efficient  and  appropri¬ 
ate  to  duplicate  the  wave  motions  that  are  governed  by 
the  time-dependent  Maxwell  equations  [1,2,3].  This  nu¬ 
merical  scheme  is  derived  from  the  eigenvalue  and  the 
eigenvector  structure  of  the  hyperbolic  equations  sys¬ 
tem.  The  procedure  is  to  diagonalize  the  three-  dimen¬ 
sional  governing  equations  into  three  one-dimensional 
Riemann  problems  [3, 4, 5, 6].  Although  this  new  numeri¬ 
cal  procedure  has  potential  to  reduce  the  required  com¬ 
puting  resources  by  allowing  larger  time  steps  and  fewer 
discretized  mesh  points  in  CEM  simulations,  substantial 
progress  in  CEM  for  practical  applications  can  finally  be 
achieved  by  incorporating  the  massively  parallel  com¬ 
puting  technique  to  deliver  the  needed  computing  re¬ 
sources. 

Recently,  through  remarkable  progress  in  mi¬ 
crochip  and  interconnect  data  link  technology,  a  host 
of  single  address,  shared  memory,  and  multiple  address 
message-passing  parallel  computers  becomes  available 
for  data  processing.  These  scalable  multi-processors 
or  multi-computers,  in  theory,  are  capable  of  providing 
essentially  unlimited  computing  resources  for  scientific 
simulations.  However,  the  effective  use  of  these  mas¬ 
sively  parallel  computers  still  rests  squarely  on  balanc¬ 
ing  the  work  load  and  keeping  the  communication  be¬ 
tween  computing  nodes  to  an  absolute  minimum  [9,10]. 
These  requirements  are  intrinsically  related  to  the  nu¬ 
merical  algorithms  and  hardware  architectures.  In  the 
present  research  effort,  attempts  were  made  to  map  a 
characteristic-based  algorithm  onto  a  message-passing 
parallel  computer  for  computational  electromagnetics. 

2  Analysis 

The  characteristic-based  fractional-step  algo¬ 
rithms  have  been  demonstrated  to  be  very  efficient  in 
solving  three-dimensional  Maxwell  equations  in  the  time 
domain  [3,4,5].  In  this  method,  the  time-dependent 
Maxwell  equations  can  be  written  in  the  flux  vector 
form  [11,12]  given  by 


d£  ^  ^  dH 

dt  ^  dx  ^  dy  ^  dz 


(1) 


^  A  partijil  differential  equations  AUxx+BUxy+CUyy-\-DU x+ 
EUy  =  F  with  —4AC  >  0,  or  all  eigenvalue  of  these  coefficient 
matrices  are  real. 


where 


The  above  partial  differential  equation  system  is 
hyperbolic,  and  constitutes  an  initial  value  problem. 
Since  the  three  coefficient  matrices  can  only  be  di¬ 
agonalized  one  at  a  time  [6,7],  the  three-dimensional 
equation  is  split  into  three  one-dimensional  formula¬ 
tions.  For  each  one-dimensional  formulation,  the  solv¬ 
ing  procedure  starts  with  solving  the  eigenvalues  and 
eigenfunctions  of  that  particular  spatial  direction.  The 
eigenvalues  for  all  three  coefficient  matrices  turn  out 
to  be  {— c,  — c,  0, 0,c,  c},  where  c  is  the  wave  propa¬ 
gating  speed.  The  signs  of  these  eigenvalues  actu¬ 
ally  control  the  direction  of  the  information  flow  as 
the  wave  propagates  through  the  computational  domain 
[3,4, 5, 6].  Based  on  the  eigenvectors,  the  field  compo¬ 
nents  are  then  expressed  in  terms  of  characteristic  vari¬ 
ables  for  that  spatial  direction.  Once  expressed  in  the 
characteristic  form,  the  one-dimensional  formulation  be¬ 
comes  that  of  the  Riemann  problem  [5,6].  After  apply¬ 
ing  the  same  procedure  to  the  other  two  dimensions, 
the  characteristic-beised  fractional-step  or  time-splitting 
method  can  now  be  summarized  in  the  following  formu- 
leis: 

10’^+'^  =  L^LyL,L,LyL^w^  (2) 


^  dWx^i  ,  ^  dWgc  i  ^ 


f  =  l,2,...,6  (3) 

f  =  l,2,...,6  (4) 

i  =  l,2,...,6  (5) 


where  WxjiVy,  andtn^  are  the  one-dimensional  character¬ 
istic  variables  associated  with  the  x,  y,  and  z  directions; 
L  is  the  one-dimensional  difference  operator;  A;  are  the 
eigenvalues  of  the  partial  differential  equation  system; 
and  n  is  the  index  of  time  level®. 


*The  fractional  Step  method  sweeps  eJl  spatial  derivatives 
twice  in  a  symmetric  sequence,  thus  advances  the  time  step 
accordingly 
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The  second-order  accurate  windward  difference 
approximation  can  be  easily  constructed  by  using  a  sin¬ 
gle  forward  or  backward  difference  approximation  ac¬ 
cording  to  the  signs  of  the  eigenvalues.  The  ability  to 
associate  the  sign  of  the  eigenvalue  with  the  direction 
of  wave  propagation  is  a  very  important  feature  of  the 
present  numerical  procedure.  The  windward  difference 
approximation  for  Lx  is  given  by 


dw 

dt 

dw 

dx 

dw 

dx 


—  111" 


2Ax 

2Ax 


-  A  <  0  (7) 
A  >  0  (8) 


where  *  takes  the  values  n  -|-  1  and  n  respectively 
for  implicit  and  explicit  procedures;  i,j,  k  are  the  indices 
of  discretized  grids.  The  windward  difference  approxi¬ 
mations  for  Ly  and  are  similar  to  that  for  Lx  ■  They 
can  be  obtained  by  replacing  x  in  Equations  (10)  and 
(11)  with  y  and  z,  and  applying  the  windward  difference 
to  the  j  and  k  indices,  respectively. 


From  Equation  (5),  the  time-splitting  solving  pro¬ 
cedure  requires  six  one  directional  numerical  sweeps  to 
update  the  solution  to  the  next  time  level.  Since  the 
formulation  for  each  numerical  sweep  is  identical,  the 
strategy  to  map  each  one  directional  sweep  to  a  mas¬ 
sively  parallel  computer  is  identical.  Thus,  a  code  based 
on  the  time-splitting  algorithm  can  be  easily  tuned  for 
maximum  performance  by  counting  and  adjusting  arith¬ 
metic  operations  for  each  numerical  sweep. 

3  Data  Partition  Schemes  and  Fine  Tuning  for 
Maximum  Parallel  Performance 

The  partition  of  data  structure  plays  a  key  role  in 
achieving  high  parallel  efficiency.  On  a  message  pass¬ 
ing  or  distributed  memory  multicomputer  system,  the 
performance  of  concurrent  computing  is  closely  tied  to 
memory  bandwidth  and  memory  latency.  Thus,  the  ba¬ 
sic  strategy  in  achieving  high  parallel  efficiency  of  any 
numerical  procedure  is  to  minimize  data  movements  be¬ 
tween  nodes.  In  this  effort,  a  computer  code  based  on 
the  characteristic-based  fractional-step  algorithm  was 
mapped  onto  the  Intel  Touchstone  Delta  at  Caltech  us¬ 
ing  different  data  partition  schemes  to  explore  the  par¬ 
allel  performance  of  the  code. 


J-D  BLOCKED  PARAUELIZATION 


Figure  1:  Hierarchy  of  data  partition 
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The  Delta  multicomputer  consists  of  a  total  of 
576  heterogeneous  nodes  [13].  They  are  designated  as 
numeric,  service,  gateway,  and  disk  nodes  to  perform 
the  computation,  frame  buffer,  network  link,  and  disk 
string  functions  respectively.  A  total  of  512  numeric 
nodes  are  available  for  parallel  computing.  Each  of 
the  numeric  nodes  has  a  peak  rate  of  80  MFlops  for 
single-precision  and  60  MFlops  for  double-precision,  and 
is  interconnected  by  a  two-dimensional  mesh  network. 
On  this  ensemble  nodal  architecture,  the  data  flow  and 
data  management  lead  to  four  different  approaches  for 
the  controlling  of  data  movements  between  nodes.  The 
most  elementary  approach  to  data  partition  is  the  one¬ 
dimensional  parallelization  in  which  the  outermost  do 
loop  of  the  numerical  sweeps  is  assigned  to  a  number 
of  numeric  nodes.  The  other  data  partition  schemes  in¬ 
clude  the  two-dimensional  page  structure  by  partition¬ 
ing  three  dimensional  space  into  cross-sectional  planes, 
the  pencil  grouping,  and  finally,  the  three-dimensional 
block  parallelization.  A  graphic  depiction  of  those  data 
partition  schemes  is  given  in  Figure  1.  The  computa¬ 
tional  domain  is  assumed  to  consist  of  IL,  JL,  and  KL 
discretized  mesh  points  in  each  direction  of  the  three 
dimensional  space;  the  computation  is  assigned  to  NP 
numeric  nodes. 

As  mentioned  in  the  beginning  of  this  section,  min¬ 
imizing  the  number  of  message  passings  is  by  far  most 
important  in  achieving  high  parallel  efficiency.  In  addi¬ 
tion,  load  balancing  is  also  very  important  in  reducing 
the  performance  disparity  among  nodes  since  the  per¬ 
formance  of  the  slowest  node  actually  determines  the 
completion  of  a  numerical  simulation.  Last,  but  not 
least,  an  ultimate  data  assigning  sequence  which  takes 
advantage  of  the  nearest  neighbor  message-passing  pri¬ 
ority  is  essential  in  reducing  the  unnecessary  message 
routing  length  and,  therefore,  improve  the  over  all  par¬ 
allel  performance.  In  the  following  discussion,  the  num¬ 
ber  of  message  passings,  load  balancing,  and  sequence 
assigning  issues  will  be  addressed. 

Although  four  data  partition  schemes  are  identi¬ 
fied  in  mapping  the  characteristic-based  code  onto  the 
Delta,  only  the  one-dimensional  parallelization  and  the 
pencil  grouping  scheme  were  investigated  in  this  effort. 
The  one-dimensional  scheme  is  the  most  straightforward 
data  partition  for  the  time-splitting,  windward  proce¬ 
dure  [5].  In  this  scheme,  the  computational  domain  of 
ILxJLxKL  mesh  points  is  divided  into  either  IL,  JL, 
or  KL  grid  planes.  Without  loss  of  generality,  the  com¬ 
putational  domain  is  assumed  to  be  divided  into  IL  grid 
planes.  Thus,  each  grid  plane  contains  JL  x  KL  mesh 
points.  The  computation  associated  with  the  JL  x  KL 
mesh  points  of  each  grid  plane  is  assigned  to  a  numeric 
node.  At  each  time  level,  each  node  performs  eight  mes¬ 


sage  passings,  which  include  messages  sent  to  and  re¬ 
ceived  from  its  nearest  four  neighboring  nodes.  In  all, 
the  total  number  of  message  passings  is  given  by  the 
value  of  2  X  (4  X  /L  —  6).  With  four  bytes  per  mes¬ 
sage,  the  message  length  associated  with  each  message 
passing  is  determined  hy  J L  x  KL  x  A.  For  example,  a 
{node  X  48  X  48)  mesh  system  will  have  a  message  length 
of  9,216  bytes. 

Since  each  node  performs  the  same  amount  of 
arithmetic  operations  in  the  one-dimensional  partition, 
load  balance  is  automatically  achieved.  This  type  of 
data  partition  can  be  viewed  as  a  form  of  task  parti¬ 
tioning. 

In  one-dimensional  parallelization  scheme,  the 
number  of  numeric  nodes  required  is  equal  to  the  num¬ 
ber  of  grid  planes,  ’IL’.  To  locate  the  desired  number 
of  numeric  nodes,  a  Delta  partition  {m,  n)  with  m  rows 
and  n  columns  should  contain  ’IL’  nodes.  Every  node  in 
this  partition  has  a  logic  number  assigned  to  it  which  a 
programmer  can  access  through  the  function  call  ”myn- 
ode()”.  The  returned  value  of  ”mynode()”  ranges  from 
0  to  nm  —  1  for  the  (m,  n)  partition.  Since  IL  equals 
nm  —  1,  an  easy  way  of  mapping  the  grid  planes  to  the 
nodes  is  to  assign  grid  plane  number  one  to  node  0,  grid 
plane  number  two  to  node  1,  and  etc.  In  other  word,  for 
the  ith  grid  plane,  the  computation  associated  with  that 
grid  plane  is  performed  by  numeric  node  i  —  1.  Unfor¬ 
tunately,  this  logical  sequence  assignment  does  not  take 
advantage  of  the  nearest  neighbor  message-passing  pri¬ 
ority  and  incurs  unnecessary  delay  in  data  movement. 
Particularly,  the  message  passing  in  Delta  is  row  biased. 
If  a  grid  plane  is  located  at  or  next  to  a  boundary  of 
the  Delta  partition,  the  messages  of  this  grid  plane  will 
have  to  be  routed  through  ‘n’  or  ‘n  —  V  nodes  before 
reaching  its  destination  grid  plane  located  in  a  different 
row  of  the  Delta  partition. 

One  way  to  avoid  this  delay  is  to  use  the  serpentine 
sequence  in  assigning  the  grid  planes  to  the  numeric 
nodes,  such  that  any  message  need  not  travel  more  than 
two  nodes  to  reach  its  destination  grid  plane  [11].  The 
serpentine  arrangement  can  be  easily  achieved  by  the 
following  program  instruction: 

.  _  J  mynode()  -f  1  x  =  even  number 

(2*a;-|-l)!t!n  —  mynode()  x  =  odd  number 

where  the  row  number  x  =  0,1,2, . m—  1.  The  ex¬ 

ample  of  a  (4, 5)  Delta  partition  with  its  physical  layout 
and  logical  node  numbers  given  by  mynode()  is  shown 
in  Table  1  (a).  The  indices  of  the  grid  planes  mapped  to 
the  (4, 5)  partition  using  logical  sequence  and  serpentine 
sequence  are  shown  in  table  1  (b)  and  (c),  respectively. 

In  the  pencil  grouping  scheme  employed  here,  the 
calculations  associated  with  two  dimensions  of  the  com- 
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Table  1:  (a)  Logical  node  numbers  of  a  (4,5)  Delta  partition  given  by  mynode().  (b)  The  grid  planes  mapped  to  the 
(4, 5)  partition  using  Logical  sequence,  (c)  The  grid  planes  mapped  to  the  (4, 5)  partition  using  Serpentine  sequence. 


putational  domain,  for  instance  IL  x  JL,  is  assigned 
to  an  {m,  n)  partition  on  the  Delta.  In  other  words, 
each  numeric  node  performs  the  calculations  associated 
with  {IL  X  JL  X  KL)/mn  grid  points.  If  the  value 
of  {IL  X  JL  X  KL)/mn  turns  out  to  be  an  integer, 
each  node  will  perform  the  same  amount  of  arithmetical 
operations  and  load  balance  is  automatically  achieved. 
If  the  value  is  non-integer,  some  nodes  will  have  more 
grid  points  to  calculate  than  some  others.  In  this  case, 
an  optimum  assignment  must  be  sought  to  balance  the 
load.  Although  in  some  cases  load  balance  is  not  auto¬ 
matically  achieved,  a  pencil  grouping  scheme  provides 
greater  flexibility  than  the  one-dimensional  scheme  in 
assigning  grid  points  to  the  numeric  nodes. 

Since  in  the  pencil  grouping  scheme,  only  the  data 
near  the  surface  of  the  pencil  must  be  transferred  to  the 
adjacent  pencils,  the  data  flow  is  reduced  compared  with 
the  one-dimensional  partition  scheme.  Thus,  a  larger 
pencil  size  will  reduce  the  overall  data  flow  and  enhance 
the  concurrent  performance.  The  number  of  message 
passings  depends  on  the  size  of  each  pencil  partition. 

4  Scope  of  Numerical  Experiments 

All  numerical  solutions  reported  in  this  paper  were 
processed  on  the  Delta  system.  The  numerical  simu¬ 
lations  were  focused  on  a  three-dimensional  oscillating 
electric  dipole.  Although  the  physics  of  the  oscillating 
dipole  is  well  known,  it  is  chosen  because  its  theoretical 
solution  [11,12]  can  be  used  for  validation.  Besides,  the 
exact  solution  contains  a  singular  behavior  at  the  cen¬ 
ter  of  the  dipole.  The  leading  term  of  this  singularity 
appears  as  the  inverse  cubic  power  of  the  radial  dis¬ 
tance  from  the  dipole  [12].  Therefore,  it  offers  a  serious 
challenge  to  the  accuracy  of  the  numerical  simulation 
and  the  robustness  of  the  solving  procedure.  To  deter¬ 
mine  the  suitable  mesh  spacing  to  capture  the  singular 
behavior  of  the  dipole,  a  three-level  mesh  spacing  refine¬ 
ment  study  was  conducted.  The  numerical  results  were 
compared  with  the  theoretical  solutions  to  assess  the  ac¬ 
curacy  of  the  simulations.  Details  comparisons  will  be 
reported  in  the  next  section. 

In  addition  to  the  validation  with  the  exact  solu¬ 


tions,  the  results  generated  by  the  parallel  codes  for  the 
one-dimensional  parallelization  and  the  pencil  grouping 
schemes  on  a  (48  x  48  x  48)  mesh  system  were  compared 
with  the  previous  numerical  results  generated  on  a  single 
vector  processor  [5].  The  parallel  results  are  identical  to 
the  serial  results. 

As  mentioned  before,  the  three  keys  to  improve 
parallel  performance  is  to  balance  the  load,  shorten  mes¬ 
sage  routine  path  lengths,  and  minimize  the  number  of 
message  passings  between  nodes.  Since  for  ID  paral¬ 
lelization,  the  load  balance  is  already  achieved,  perfor¬ 
mance  evaluation  was  concentrated  on  the  latter  two. 
In  the  next  section,  parallel  performance  of  the  logical 
and  serpentine  sequence  mapping  schemes  will  be  com¬ 
pared  using  mesh  systems  of  {node  x  48  x  48),  where 
node  equals  4,  8,  16,  32,  64,  256,  384,  and  512;  per¬ 
formance  comparison  for  different  number  of  message 
passings  will  be  carried  out  using  {node  x  48  x  48)  and 
{node  X  96  X  96)  mesh  systems.  For  the  pencil  group¬ 
ing  parallelization  scheme,  the  three  keys  to  improve 
parallel  performance  still  apply.  The  pencil  grouping 
performance  will  be  reported  as  well. 

All  the  parallel  performances  were  evaluated  based 
on  the  timing  data  obtained  from  the  execution  phase 
and  data  output  phase  of  the  codes.  The  execution 
phase  includes  updating  the  fields  at  all  grid  points  from 
the  first  time  level  up  to  the  last  time  level  of  the  simu¬ 
lation;  the  data  output  phase  includes  writing  the  calcu¬ 
lated  field  results  and  the  coordinates  of  the  grid  system 
to  the  disk  for  post  processing.  The  program  initializa¬ 
tion  time  of  the  present  code  is  negligible  compared  to 
the  program  execution  and  data  output  time;  therefore, 
it  will  not  be  reported  here.  (In  the  following  text,  data 
output  is  sometime  referred  to  as  data  I/O  since  the 
output  performance  is  tied  into  the  I/O  performance  of 
the  computer  system.) 

The  execution  times  were  collected  at  the  480th 
time  level.  The  duration  of  480  time  steps  was  selected 
for  all  calculations  such  that  the  wave  traversed  more 
than  ten  times  the  distance  across  the  entire  computa¬ 
tional  domain.  This  time  duration  is  deemed  sufficient 
for  the  numerical  procedure  to  yield  meaningful  results 
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and  repeatable  timing  data.  From  the  execution  time, 
the  data  processing  rate  and  scale  factor  are  derived. 
At  each  time  level  and  grid  point,  both  one-dimensional 
and  the  pencil  grouping  approach  must  perform  at  least 
492  mixed  integer  and  real  number  arithmetic  opera¬ 
tions.  The  data  processing  rate  is  computed  by  the  total 
number  of  arithmetic  operations  performed  during  the 
480  time  steps  divided  by  the  maximum  and  the  mini¬ 
mum  time  elapsed  of  all  nodes  in  use.  These  minimum 
and  maximum  values  bracket  the  range  of  performance 
among  nodes.  In  other  word,  the  data  processing  rate 
in  Gigaop  (10®  arithmetic  operations  per  second)  is  cal¬ 
culated  by  multiplying  the  numbers  of  iterations,  arith¬ 
metic  operations,  and  total  grid  points,  then  dividing 
by  the  executing  time  (480  x  492  x  (IL  x  JL  x  KL)/ 
(execution  time  xlO®)).  This  calculated  data  rate  is  a 
conservative  estimate  becanse  the  number  492  does  not 
include  47  basic  library  function  calls  and  nine  system 
message  passing  calls  at  every  time  step.  These  library 
function  calls  include  mostly  trigonometry  and  square 
root  calculations.  Together  with  the  system  message 
passing  calls,  they  may  account  for  a  significant  amount 
of  the  execution  time.  The  scale  factor  of  the  code  is 
defined  by  the  time  elapsed  for  each  array  of  nodes  used 
then  normalized  by  the  execution  time  required  for  four 
nodes.  The  execution  time  for  the  four-node  calculation 
was  used  as  timing  basis  because  in  the  present  code,  at 
least  four  grid  points  are  required  to  specify  the  bound¬ 
ary  conditions  at  the  dipole  and  at  the  far  field. 

The  data  output  phase  of  the  one-dimensional 
parallelization  scheme  requires  78  arithmetic  operations 
and  38  system  message  passing  calls.  Most  of  the  38  sys¬ 
tem  calls  are  either  synchronous  or  cisynchronous  mes¬ 
sage  passing  instructions  to  generate  one  unformatted 
field  data  file  and  three  unformatted  grid  files.  The 
field  data  file  contains  all  six  components  of  the  elec¬ 
tromagnetic  field  at  every  grid  point;  and  the  grid  files 
contain  the  three  coordinates  of  the  grid  system.  The 
longest  and  the  shortest  time  periods  required  to  out¬ 
put  those  unformatted  files  for  4-node  up  to  512-node 
simulations  were  recorded.  From  these  timing  data,  the 
scale  factor  for  each  simulation  was  again  normalized  by 
the  4-node  result  to  determine  the  possible  performance 
degradation  of  data  management. 

5  Discussion  of  Results 

5.1  Mesh  refinement  study 

The  originally  planned  mesh  refinement  included 
mesh  systems  of  (48  x  48  x  48),  (96  x  96  x  96),  and  a 
finest  mesh  system  of  (192  x  192  x  192).  However,  as  the 
number  of  processors  increased  beyond  96,  the  disparity 
in  execution  times  among  processors  grew  to  a  factor  of 
11.6  for  108  nodes,  and  14.43  for  144  nodes.  In  other 


words,  the  observed  scalable  performance  of  the  present 
computer  program  for  the  cases  of  (48  x  48  x  48)  and 
(96  x  96  X  96)  was  not  sustainable  for  the  mesh  system 
of  (192  X  192  X  192).  Therefore,  the  mesh  refinement 
study  was  forced  to  reduce  the  scope. 

Three  mesh  systems  of  (24  x  24  x  24) ,  (48  x  48  x  48) , 
and  (96  x  96  x  96)  were  used  for  the  numerical  resolu¬ 
tion  study  instead.  For  this  reduced  scope  of  numerical 
experiment,  the  coarse  mesh  system  has  only  24  grid 
points  in  each  coordinate.  The  mesh  point  density  is 
sufficient  to  resolve  the  wave  motion  [1,2]  away  from  the 
dipole,  but  it  may  be  deficient  in  simulating  the  singular 
field  behavior  near  the  dipole  [11,12]. 

The  mesh  refinement  results  for  the  radial  com¬ 
ponent  of  the  electric  field  are  presented  in  Figure  2. 
The  calculation  from  the  (24  x  24  x  24)  mesh  system 
revealed  significant  numerical  oscillations  in  an  attempt 
to  overcome  the  large  truncation  error  near  the  dipole. 
As  shown  in  Figure  2,  the  results  exhibits  a  steady  im¬ 
provement  as  the  mesh  systems  becomes  finer.  How¬ 
ever,  the  finest  mesh  is  still  not  fine  enough  to  overcome 
the  stringent  demand  for  accurate  simulation  near  the 
dipole.  Nevertheless,  the  present  numerical  procedure 
has  demonstrated  the  robustness  in  treating  the  prob¬ 
lem  containing  a  singular  behavior. 
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Figure  2:  Comparison  of  computed  radial  components 
of  the  electric  fields  with  theory 

Figure  3  depicts  the  numerical  results  of  the  cir¬ 
cumferential  component  of  electric  intensity  on  the  three 
mesh  systems.  The  numerical  result  obtained  on  the 
finest  mesh  attained  the  best  agreement  with  the  ana¬ 
lytical  result.  The  maximum  deviation  from  the  theoret¬ 
ical  results  is  merely  a  fraction  of  one  percent.  Another 
desired  feature  of  the  characteristic-based  formulation 
also  stands  out.  At  a  non-dimensionalized  radial  dis- 
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tance  of  0.14  and  beyond,  all  three  solutions  agree  well 
with  the  theory  and  show  no  wave  reflection  from  the 
truncated  numerical  boundary. 


r 

Figure  3:  Comparison  of  computed  circumferential  com¬ 
ponents  of  the  electric  flelds  with  theory 

The  comparison  of  magnetic  azimuthal  component 
of  the  three  mesh  systems  is  given  in  Figure  4.  Although 
the  leading  term  singularity  of  the  magnetic  field  has  a 
lower  order  asymptote  than  that  of  the  electric  field,  the 
numerical  solutions  of  the  magnetic  field  generated  on 
the  three  mesh  systems  behave  similarly  to  those  of  the 
circumferential  electric  field  component.  The  numerical 
resolution  produced  by  the  finest  mesh  system,  however, 
is  able  to  suppress  the  numerical  oscillation  near  the 
dipole.  Finally,  the  reflected  wave  from  the  truncated 
numerical  boundary  is  completely  absent. 


Figure  4:  Comparison  of  computed  azimuthal  compo¬ 
nents  of  the  magnetic  fields  with  theory 


5.2  Serpentine  sequence  vs.  logical  sequence 


In  this  experiment,  the  one-dimensional  paral¬ 
lelization  scheme  and  the  grid  system  {node  x  48  x  48) 
were  used  so  that  every  node  performs  exactly  the 
same  amount  of  computations  and  writes  out  the  same 
amount  of  outputs. 

Table  2  shows  the  data  processing  rates  in  Gi- 
gaops.  The  performance  difference  between  the  logical 
sequence  and  the  serpentine  message  path  arrangements 
begins  to  appear  for  the  number  of  nodes  beyond  128. 
Since  the  serpentine  arrangement  took  advantage  of  the 
nearest  neighbor  priority  in  message  passing  hierarchy, 
the  performance  degradation  is  less  than  that  of  the  log¬ 
ical  sequence.  In  fact,  to  complete  a  512-node  compu¬ 
tation  using  serpentine  procedure,  the  slowest  node  re¬ 
quired  only  3.1  percent  more  execution  time  than  the 
fastest  node;  whereas,  in  the  logical  sequence  arrange¬ 
ment,  the  slowest  node  is  10.4  percent  slower  compared 
with  the  fastest  node.  For  the  same  512-node  simula¬ 
tion,  the  fastest  node  operating  on  the  serpentine  ar¬ 
rangement  attained  a  data  processing  rate  of  5.966  Gi- 
gaops;  while  the  fastest  node  on  the  logical  sequence 
arrangement  only  attained  5.412  Gigaops. 
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Figure  5:  Data  processing  rates  on  the  Delta 


The  data  processing  rates  of  the  serpentine  and 
logical  sequences  are  plotted  in  Figure  5.  The  improved 
data  processing  rate  of  the  serpentine  sequence  over  that 
of  the  logical  sequence  is  clearly  demonstrated.  As  can 
be  seen  from  the  figure,  the  performance  of  the  logical 
sequence,  which  did  not  honor  the  data  passing  priority 
of  the  nearest  neighbor  hierarchy,  started  to  lag  behind 
that  of  the  serpentine  sequence  when  more  than  128 
numeric  nodes  were  used.  From  this  experiment,  the 
serpentine  procedure  appears  to  be  more  effective  for 
the  distributed  memory  system.  However,  for  both  se¬ 
quences,  an  increasing  disparity  of  the  data  processing 
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Number  of  nodes 

Serpentine 

Maximum  Minimum 

Logical  Sequence 
Maximum  Minimum 

4 

0.046 

0.046 

0.046 

0.046 

8 

0.092 

0.092 

0.092 

0.092 

16 

0.184 

0.184 

0.184 

0.169 

32 

0.368 

0.368 

0.369 

0.369 

64 

0.738 

0.738 

0.736 

0.707 

128 

1.472 

1.472 

1.445 

1.430 

246 

2.967 

2.952 

2.829 

2.798 

384 

4.490 

4.413 

4.105 

4.028 

512 

5.966 

5.704 

5.412 

5.274 

Table  2:  Comparison  of  data  rates  between  serpentine  and  logical  sequences  of  an  (node  x  48  x  48)  mesh  system 


rate  among  nodes  was  noted  as  the  number  of  nodes 
increased  beyond  64. 

The  scale  factors  of  the  two  mapping  schemes  are 
demonstrated  in  Figure  6.  All  data  processing  rates  are 
normalized  by  the  value  of  the  4-node  simulation.  The 
superior  scalable  property  of  the  serpentine  sequence 
over  that  of  the  logical  sequence  is  obvious.  The  scala¬ 
bility  of  the  present  code  up  to  512  nodes  is  perceived 
within  a  performance  degradation  of  less  than  3.1  per¬ 
cent.  In  fact,  the  significant  degradation  only  appeared 
when  all  512  numeric  nodes  were  employed. 
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Figure  6:  Scalability  of  executing  time  on  the  Delta 

Figure  7  presents  the  data  output  time  in  seconds. 
The  data  output  time  varied  widely  from  0.24  to  32.20 
seconds  for  the  serpentine  sequence,  and  from  0.28  to 
162.80  seconds  for  the  logical  sequence.  The  shortest 
waiting  time  over  the  entire  range  of  nodes  used  was 
only  0.24  seconds.  Unfortunately,  the  aggregated  time 
required  by  the  slowest  node  in  computation  and  data 
output  actually  determines  the  completion  of  a  numeri¬ 
cal  simulation.  From  Figure  7,  the  data  output  time  in¬ 


creases  almost  linearly  with  the  number  of  nodes  in  use. 
For  the  512-node  calculations,  the  serpentine  sequence 
used  32.2  seconds  while  the  logical  sequence  needed 
162.8  seconds  to  output  the  same  amount  of  data. 
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Figure  7:  Data  I/O  time  on  the  Delta 

Figure  8  depicts  the  widely  varying  data  I/O  per¬ 
formance  among  nodes  for  the  entire  range  of  available 
nodes.  Again  the  performance  degrades  rapidly  as  the 
number  of  nodes  in  use  increases.  The  logical  sequence 
yields  the  maximum  parallel  performance  discrepancy 
among  nodes.  The  ratio  between  the  most  and  least 
efficient  nodes  is  as  high  as  580.4.  The  serpentine  se¬ 
quence  reduces  the  disparity  to  a  value  about  134.2.  If 
this  behavior  is  a  common  trend  for  all  distributed  mem¬ 
ory  computer  systems,  this  deficiency  shall  be  a  pacing 
item  for  research  in  concurrent  computing. 

Concurrent  performance,  similar  to  those  dis¬ 
cussed  above,  were  also  observed  for  (96  x  96  x  96)  and 
(96  X  192  X  192)  computations.  The  data  processing 
rate  per  node  was  maintained  in  the  range  of  12.89  to 
11.65  Megaops.  These  data  processing  rates  are  far  be- 
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low  the  nominal  60  to  80  Megaops  performance  level  of 
the  i860  microprocessor  [13].  Additional  performance 
improvement  of  the  present  numerical  procedure  using 
the  Delta  system  is  still  possible. 
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Figure  8;  Non-scalability  of  data  I/O  on  the  Delta 


5.3  Comparison  between  {node  x  48  x  48)  and 

{node  X  96  X  96)  mesh  systems 

For  the  one-dimensional  parallelization  scheme, 
the  vector  length  is  doubled  and  the  message  passing 
length  is  quadrupled  from  the  coarse  grid  of  {node  x 
48  X  48)  to  the  refined  grid  of  {node  x  96  x  96).  Since 
the  serpentine  mapping  scheme  was  proven  to  give  bet¬ 
ter  parallel  performance,  it  was  used  for  the  experiment 
reported  in  this  subsection. 

For  the  {node  x  48  x  48)  systems,  the  execution 
time  for  each  node  varied  from  45.61  to  52.90  seconds 
regardless  of  whether  the  computations  were  conducted 
on  4  or  512  numeric  nodes.  For  the  {node  x  96  x  96) 
case,  the  execution  time  ranges  from  158.37  to  173.16 
seconds.  The  data  processing  rates  based  on  the  above 
execution  time  is  given  in  Table  3.  The  maximum  data 
processing  rate  of  the  refined  mesh  system  is  9.8  per¬ 
cent  greater  than  the  coarser  mesh  computation;  and 
the  highest  achieved  value  was  6.551  Gigaops.  Appar¬ 
ently  at  a  message  passing  length  of  36,864  bytes  and 
8  message  passings  per  node  per  time  step,  the  higher 
data  processing  rate  by  a  greater  vector  length  has  not 
been  hampered  by  the  limiting  memory  bandwidth.  As 
a  result,  the  deviation  between  the  fastest  and  the  slow¬ 
est  execution  rates  is  less  than  one  percent. 

The  execution  times  expressed  in  scale  factors  for 
the  {node  x  96  x  96)  and  {node  x  48  x  48)  systems  are 
shown  in  Figure  9.  Small  variations  of  data  process¬ 
ing  rates  for  both  grid  systems  were  observed  over  the 


entire  range  of  nodes  used.  The  results  of  the  coarse 
grid  computations  exhibited  a  wider  performance  de¬ 
viation  among  nodes  than  the  finer  grid  calculations. 
The  latter,  however,  also  degraded  noticeably  from  the 
4-node  to  the  8-node  result.  The  difference  is  about 
4.6  percent.  As  the  number  of  nodes  increases  to  512, 
the  degradation  of  scalability  of  the  present  procedure 
reached  the  maximum  value  of  7.6  percent.  But  the 
major  portion  of  the  performance  deficiency  is  incurred 
at  the  4-node  to  8-node  transition.  In  all,  the  present 
algorithms  mapped  to  Delta  system  has  demonstrated 
an  acceptable  scalability  of  data  processing  rate  up  to  a 
grid  system  consisting  of  512  x  96  x  96  grid  points. 
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Figure  9:  Execution  time  scale  factors  for 
the  coarse  and  fine  mesh  systems 

The  data  output  time  varies  from  0.2369  to  32.20 
for  the  coarse  grid  system  and  from  0.1324  seconds  for 
the  16-node  calculation  to  30.02  seconds  for  the  512- 
node  simulation  for  the  finer  grid  system.  Figure  10 
shows  the  scale  factors  of  the  output  timing  results  for 
the  coarse  and  refined  mesh  systems.  Clearly,  the  I/O 
performance  degrades  nearly  linearly  as  the  number  of 
nodes  in  use  increases.  The  worst  I/O  performance  oc¬ 
curred  when  all  512  nodes  were  engaged.  In  this  case, 
the  fastest  and  the  slowest  nodes  used  11.51  and  226.8 
times,  respectively,  the  time  required  to  output  the  same 
amount  of  data  eis  the  fastest  node  of  the  16-node  com¬ 
putation.  This  observation  further  asserts  that  the  scal¬ 
able  I/O  performance  shall  be  a  pacing  item  for  research 
in  concurrent  computing. 

5.4  Pencil  grouping  scheme 

Owing  to  similar  peculiarities  of  the  Delta  com¬ 
piler  observed  for  the  one-dimensional  parallelization, 
only  a  limited  number  of  timing  data  from  the  pencil 
partition  is  obtained  at  present  time.  In  short,  the  ser¬ 
pentine  sequence  did  improve  the  performance  over  that 
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Number  of  nodes 

node  X  96  X  96 
Maximum  Minimum 

node  X  48  X  48 
Maximum  Minimum 

4 

0.055 

0.055 

0.046 

0.046 

8 

0.105 

0.105 

0.092 

0.092 

16 

0.207 

0.207 

0.184 

0.184 

32 

0.413 

0.413 

0.368 

0.368 

64 

0.823 

0.822 

0.738 

0.738 

128 

1.616 

1.609 

1.472 

1.472 

246 

3.305 

3.278 

2.967 

2.952 

384 

4.915 

4.884 

4.490 

4.413 

512 

6.551 

6.505 

5.966 

5.704 

Table  3:  Comparison  of  data  rates  between  node  x  96  x  96  and  node  x  48  x  48  mesh  systems 


Delta 

(m,n) 

192x1 

96x2 

64x3 

48x4 

32x6 

IL  X  JL  (Pencil  Dimension) 
24x8  16x12  12x16  8x24 

6x32 

4x48 

3x64 

2  x96 

1x192 

(6,32) 

13.4 

9.73 

10.6 

10.7 

12.1 

13.1 

14.2 

15.7 

16.8 

16.8 

16.7 

14.3 

13.2 

8.94 

(8,24) 

12.6 

10.6 

11.6 

12.1 

13.0 

13.9 

15.7 

16.8 

17.3 

16.4 

16.9 

14.0 

13.2 

9.55 

(12,16) 

12.5 

10.6 

12.0 

13.1 

14.0 

15.7 

17.0 

18.2 

17.3 

16.8 

16.5 

14.2 

13.3 

9.75 

(16,12) 

13.1 

11.3 

13.1 

13.9 

15.4 

16.2 

18.4 

17.5 

17.1 

16.2 

14.3 

13.2 

13.2 

10.1 

Table  4:  Execution  rates  (Megaops)  for  a  (192  x  192  x  192)  grid  system  by  pencil  grouping  scheme 


of  the  logical  sequence  as  the  one-dimensional  case.  The 
timing  information  for  a  (192  x  192  x  192)  mesh  system 
mapped  to  the  Delta  using  pencil  grouping  scheme  with 
various  combinations  of  pencil  dimension  (ILxJL)  and 
Delta  partition  (m,  n)  is  tabulated  in  Table  4.  The  vari¬ 
ation  of  data  processing  rates  among  the  different  pencil 
partitions  for  a  fixed  Delta  partition  is  significant.  In 
general,  slower  performance  is  observed  when  the  pencil 
dimension  and  Delta  partition  are  at  the  extreme  sit¬ 
uations  when  the  pencil  grouping  scheme  degenerates 
into  one- dimensional  parallelization.  For  example,  the 
pencil  dimensions  of  (192  x  1)  or  (1  x  192)  on  a  Delta 
partition  of  (6,32)  has  a  slower  data  processing  rate  of 
13.2  or  8.94  Megaops,  respectively. 

A  general  performance  improvement  for  an  {IL  x 
JL)  pencil  dimension  is  noted  when  IL  nearly  equals 
JL.  The  best  performance  for  a  given  Delta  partition  is 
observed  when  IL  equals  m  and  JL  equals  n,  as  under¬ 
lined  in  Table  4.  The  timing  data  of  the  pencil  partition 
has  revealed  a  parallel  efficiency  gain  of  37%  over  that 
of  the  one-dimensional  partition. 

6  Conclusions 

A  fractional-step,  upwind  numerical  procedure  for 
solving  the  three-dimensional,  time-domain  Maxwell 
equations  has  been  ported  to  the  Intel  Touchstone  Delta 
multicomputer.  The  concurrent  computations  dupli¬ 
cated  the  results  from  earlier  numerical  results  and 


compared  well  with  theory.  For  the  mesh  systems  of 
(node  X  48  X  48)  and  (node  x  96  x  96),  the  numerical 
procedure  is  scalable  up  to  512  numeric  nodes  with  only 
up  to  7.6  percent  performance  degradation.  The  fastest 
data  processing  rate  is  6.551  Gigaops  and  the  sustained 
overall  performance  is  clocked  at  5.704  Gigaops.  Further 
increased  data  processing  rate  is  still  possible. 
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Figure  10:  I/O  scale  factors  for  the  coarse 
and  fine  mesh  systems 

The  scalability  of  concurrent  computing  is  sus¬ 
tained  up  to  a  simulation  eight  times  the  size  of  the 
baseline  case.  The  scalable  performance  failed  at  the 
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fourth  level  of  grid  point  enrichment  (192  x  192  x  192). 
Although  the  architecture  of  the  Touchstone  Delta  mul¬ 
ticomputer  and  its  usefulness  are  impressive,  consistent 
performance  in  scaling  up  for  massive  data  bases  re¬ 
mains  as  a  necessary  research  emphasis.  The  scalable 
data  1/0  is  also  identified  as  a  pacing  item  for  intense 
research  for  attaining  high  performance  computation. 

The  one-dimension  parallelization  has  been 
shown  as  a  suitable  data  partition  procedure  for  a 
characteristic-based,  windward  difference  algorithm  in 
solving  the  time  dependent,  three  dimensional  Maxwell 
equations.  Further  parallel  efficiency  improvement  by 
pencil  structure  shows  developable  potential. 
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ABSTRACT:  The  application  of  the  Biconjugate  Gradient 
FFT  method  to  the  thin  co^ucting  plate  problem  is 
investigated.  Upon  comparing  with  a  Conjugate  Gradient 
FFT  method,  it  is  found  that  the  Biconjugate  Gradient 
solution  requires  a  relatively  larger  error  tolerance  to 
achieve  a  comparatively  well-behaved  current 
distribution.  Coupled  with  the  requirement  of  only  one 
matrix-vector  product  per  iteration,  the  computational 
cost  of  the  Biconjugate  Gradient  method  is,  thus,  much 
smaller  than  those  previously  reported  in  the  literature.  Of 
particular  importance  is  the  use  of  the  incident  electric 
field  as  a  starting  estimate  to  alleviate  the 
non-convergence  behaviour  which  is  usually  associated 
with  the  apfdication  of  a  Biconjugate  Gradient  approach  to 
conducting  plates  at  grazing  angle.  For  other  angles  of 
incidence,  it  is  shown  that  this  procedure  also  accelerates 
the  resulting  convergence  rates  as  compared  to  those 
obtained  by  simply  using  zero  as  an  initial  estimate. 

1.  INTRODUCTION 

The  application  of  FFT-based  methods  for  the  flat, 
conducting  plate  problem  has  been  receiving  a  lot  of 
attention  in  the  computational  electromagnetics 
community  since  the  late  1980s.  Most  of  these  use  ^e 
Conjugate  Gradient  method  (CGM)  [1]  as  the  iterative 
algoridim  in  conjunction  with  the  fast  Fourier  transform 
(ITT)  to  solve  for  the  discretized  matrix  equation  [2]-[6]. 
This  combination,  commonly  known  as  the  CGFFT 
method,  is  attractive  both  in  the  characteristics  of 
convergence  associated  with  the  CGM  and  the  reduced 
computational  costs  associated  with  using  the  FFTs  to 
approximate  all  occurring  convolutions. 

The  finite-termination  property  of  the  CGM  is,  however, 
only  valid  when  the  defining  matrix  (or  operator)  of  the 
linear  system  to  be  solved  is  Hermitian  and 
positive-definite  (HPD).  For  other  types  of  systems,  it  is 
necessary  to  premultiply  (implicitly)  both  sides  of  the 

pertinent  system  by  the  Hermitian  version  A*  of  the  matrix 
A  so  that  the  resulting  matrix  is  HPD.  The  latter  system, 
often  referred  to  as  the  normal  equation,  is,  however,  more 
ill-conditioned  than  the  origin^  system,  which,  in  turn, 
makes  the  convergence  rate  of  the  resulting  algorithm 
slower  than  that  associated  with  the  original  system.  For 
this  reason,  it  is  worthwhile  exploring  iterative  algorithms 
that  can  be  applied  directly  to  the  system  of  equations  in 
mind. 


One  such  method  is  the  Biconjugate  Gradient  (Bi-CGM) 
metfiod,  first  developed  by  Lanczos  for  finding  the 
eigenv^ues  of  an  unsymmetric  system  [7]  and  later 
extended  by  Fletcher  [8]  and  Jacobs  [9]  to  treat  complex 
indefinite  systems.  When  applied  to  a  symmetric  system, 
the  Bi-CGM  has  one  more  advantage  over  the  CGM  in  that 
only  one  matrix-vector  operation  is  needed  within  each 
iteration  step  [10].  The  combination  of  the  Bi-CGM  with 
the  FFT  to  solve  a  convolutional  matrix  equation  is  also 
possible  in  exactly  the  same  way  as  a  conventional  CGFFT 
formulation.  Despite  this,  there  are  only  a  few  applications 
of  the  Bi-CGFFT  method  to  electromagnetic  scattering 
problems  in  the  literature  compared  with  the  abundance  of 
solutions  obtained  by  its  CGFFT  counterpart. 
works  on  the  performance  of  the  Bi-CGFFT  fw 
conducting  problems  have,  however,  put  the  Bi-CGFFT 
in  a  very  favourable  position  compared  with  a  conventional 
CGFFT  approach  [11]-[12]. 

In  this  paper,  the  paformance  of  the  Biconjugate 
Gradient  FFT  method  when  applied  to  a  thin  conducting 
plate  under  plane-wave  incidence  is  investigated.  The 
non-convergence  behaviour  which  usually  happens  when 
the  Bi^GFFT  is  applied  to  a  conducting  plate  at  grazing 
incidence  is  overcome  by  choosing  the  initial  guess  to  be 

the  incident  electric  field  E‘.  A  systematic  study  is  also 
conducted  of  the  dependence  of  the  Bi-CGFFT  method  on 
the  use  of  this  initial  estimate  procedure  for  other  angles  of 
incidence.  Greatly  accelerated  rates  of  convergence  is 

shown  to  result  when  E‘  is  used  as  an  estimate  for  the 
initial  unknown  current  distribution  instead  of  simply 
using  zero  as  commonly  adopted  in  most  FFT-based 
implementations.  The  layout  of  the  paper  is  as  follows. 
Section  2  shows  the  functional  form  of  the  Bi-CGM  when 
ajqrlied  to  a  symmetric,  indefinite  system.  This  is  followed 
by  Section  3  which  deals  with  the  numerical  solutions 
obtained  by  the  Bi-CGFFT  and  the  CGFFT  methods  for 
a  variety  of  incident  configurations.  Finally,  Section  4 
draws  the  conclusion  on  the  applicability  of  the  Bi-CGFFT 
based  on  the  results  obtained  in  Section  3. 

2.  THE  BICONJUGATE  GRADIENT  METHOD 
FOR  SYMMETRIC  COMPLEX  SYSTEMS 

Table  1  shows  the  computationally  compact  form  of  the 
Biconjugate  Gradient  method  [10]  when  applied  to  a 
symmetric  system  of  the  form 
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AJ  =  E‘.  (1) 

where  J  denotes  the  unknown  current  distribution  to  be 
solved  for.  The  impedance  mabix  A  on  the  left-hand  side 
of  (1)  represents  the  couplings  between  the  diffi^nt  cells 
on  the  plate  and  is  often  formulated  in 
discrete-convolutional  form  to  enable  die  FFT  to  be 
efificiendy  applied  whenever  there  is  a  matrix-vector  to  be 
computed  [6]. 


Apart  from  the  definition  of  the  inner  product  (to 
calculate  and  the  form  of  the  Biconjugate 

Gradient  method  presented  in  Table  1  is  almost  identical  to 
a  classical  Conjugate  Gradient  formulation  for  a  real  and 
symmetric  positive-definite  (SPD)  system  [1].  In  both 
methods,  only  one  matrix-vector  product  in  the  form  of 

Ap(")  is  required  at  each  iteration  step.  For  complex 
indefinite  systems,  however,  the  CGM  has  to  be  modified 
to  avoid  the  possible  division  by  zero  and  substantial  errcw 
growth  when  this  situation  is  nearly  reached  [8].  The  main 
resulting  cost  is  an  extra  matrix-vector  product  in  the  form 

of  where  A*  denotes  the  complex  conjugate  of  the 
matrix  A  [1],  [10].  Thus,  the  computational  cost 
^sociated  with  a  symmetric  indefinite  Bi-CGM  ^proach 
is  expected  to  be  half  of  those  required  in  a  CGM 
formulation  when  ^plied  to  the  same  problem^. 

The  convergence  criteria  for  a  Biconjugate  Gradient 
FFT  method  is  specified  in  the  same  way  as  for  its 

Conjugate  Gradient  counterpart.  Let  denotes  the 


1 .  For  non-symmetric  systems,  the  application  of  aBiconjugate 
Gradient  method  still  requires  two  matrix-vector  products  per 
iteration  step.  However,  as  most  EM  scattering  formulations 
are  qmimetric  in  nature,  the  algorithm  with  r^uced  costs  is 
applicable. 


residual  field  at  the  n'*  step,  then  the  stopping  criteria  is 
defined  in  tmns  of  the  normalized  residual  norm, 


where  F^®)  denotes  the  initial  field  residual  and  ll  •  || 
denotes  the  norm  of  a  complex  vector.  At  each  step,  (2)  is 

COTiparedtoane/rortoterance  e  which  is  specified  by  the 
user.  The  usual  stopping  criteria  for  a  well-behaved 
CGFFT  solution  is  0.01  [5]-[6],  but  as  will  be 
demonstrated  in  the  next  section,  the  corresponding  value 
for  aBi-CGFFT  solution  can  be  much  larger,  thus  reducing 
consid^bly  the  required  number  of  iterations,  and  thus 
computational  costs. 

An  mpcnlant  paraineter  which  is  often  neglected  in 
other  Biconjugate  Gradient  formulations  is  the  form  of  the 

initial  guess  It  is  often  assumed  that  this  is  taken  to 
be  zero.  In  Section  3,  we  demonstrate  that  choosing  the 
initial  guess  to  be  the  incident  electric  field  E‘,  not  only 
alleviates  the  non-convergence  problem  that  is  usually 
^sociated  with  the  application  of  the  Bi-CGFFT  at  grazing 
incidence,  but  also  helps  to  accelerate  the  convergence 
rates  associated  with  other  angles  of  incidence. 

3.  A  COMPARATIVE  STUDY  OF 
COMPUTATIONAL  EFFICIENCIES  OF  THE 
BI-CGFFT  AND  CGFFT  METHODS  FOR  THE 
THIN  CONDUCTING  PLATE 

In  this  section,  the  performance  of  the  Bi-CGFFT 
method  when  applied  to  the  pulse-basis  formulation 
^raitly  proposed  by  Tran  and  McCowen  [6]  is 
investigate  For  each  incident  configuration,  a 
comparison  is  made  with  a  corresponding  CGFFT 
formulation.  The  fact  that  the  Bi-CGFFT  generally 
requires  half  the  computational  workload  compared  to  its 
CGFFT  counterpart  when  applied  to  a  symmetric  system  is 
well-known  and  has  been  investigated  elsewhere  [11], 
[12].  The  aim  of  this  section  is  to  show  that  the  efficiency 
of  the  Bi-^GFFT  can  be  enhanced  further  by  virtue  of  the 
fact  that  it  needs  a  relatively  larger  tolerance  to  achieve  a 
well-behaved  current  distribution  as  that  generated  by  the 
CGFFT.  The  importance  of  the  use  of  the  incident  electric 
field  as  a  starting  estimate  in  a  Bi-CGFFT  formulation  will 
also  be  investigated.  All  numerical  results  are  performed 
on  a  VAX  8820  computer. 

A.  Broadside  incidence 

Fig.  1  shows  the  broadside  current  distributions  obtained 
by  the  two  FFl  -based  methods  on  a  Uo*  Uo  square  plate 

using  mesh  sizes  Ax  =  Ay  =  0.0303Ao  with  different 
error  toletmices.  It  can  be  seen  from  this  figure  that  the 
Bi-CGFFT  requires  a  very  coarse  tolerance  to  yield  a 
well-behaved  solution,  as  compared  with  a  much  finer 
tolerance  as  required  by  the  CGFFT  to  achieve  essentially 
the  same  results.  When  applied  with  tiie  same  tolerance  as 

the  Bi-CGFFT,  i.e.  €  =0.1,  the  CGFFT  results  are 
hardly  recognizable.  For  the  x  -  component,  only  a  vague 
resemblance  of  the  expected  lobes  along  the  front  and  back 
edges  can  be  assumed  (see  [5]-[6]  for  a  discussion  of  these 
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results  with  respect  to  different  discretization  procedures); 
the  ripples  in  the  centre  of  the  scattering  plate  are  too  high 
in  comparison  with  the  edge  lobes.  The  behaviour  of  the 
y  -  componentat  the  sameerror  tolerance  is  even  moreout 
of  shape:  the  lobes  along  the  four  edges  are  nearly  of  the 
same  amplitudes,  as  are  the  ripples  in  the  centre  of  the 
scattering  plate.  Although  these  anomalies  are 

self-corrected  at  the  finer  tolerance  of  6  =  0.1  *  10"^,  the 
resulting  computational  cost  is  much  higher  compared  with 
the  corresponding  Bi-CGFFT:  the  CPU  times  required  by 


the  CGFFT  on  a  VAX  8820  is  08:40  minutes  (172 
iterations)  as  compared  with  only  01:17  minutes  (27 
iterations)  for  the  Bi-CGFFT — an  increase  of  676%.  Yet 
there  is  little  difference  between  the  two  numerical 
solutions:  the  only  noticeable  discrepancy  is  in  the 
co-polarized  y  -  component,  where  the  CGFFT  current 
density  is  slightly  smoother  along  the  four  edges.  The 
predicted  current  distribution  from  the  Bi-CGFFT  method 
for  the  dominant  j:- component  is  virtually 
indistinguishable  from  its  CGFFT  counterpart. 


(b)  y-component 

Fig.  1:  Comparison  of  the  current  distributions  by  the  CGFFT  and  Bi-CGFFT  methods  with  different  error  tolerances 
on  a  square  PEC  plate  of  size  Uo  ♦  Uo  •  Normal  incidence. 


B.  Grazing  incidence 

An  important  incident  confi^ation  of  the  plate 
scattering  problem  is  the  grazing  incidence  case.  When 
applied  to  a  geometric  shape  that  possesses  edges  ^d 
comers  such  as  a  rectangular  plate,  the  rapidly  changing 
behaviour  of  the  current  density  near  the  edge  and  the 
comer  diffraction  effects  makes  its  computation  a 
particularly  difficult  one  [13].  The  Bi-CGFFT  ^t  gives 
the  solution  for  the  broadside  incidence  in  Fig.  1  fails  to  give 
a  convergent  result  when  the  plate  is  subject  to  the  grazing 


angle.  This  undesirable  behaviour,  commonly  known  as 

the  stagnation  problem,  arises  when  the  constant 
becomes  close  to  machine-zero  and  causes  both  the  current 

density  and  field  residual  to  be 

non-incremental  as  the  iteration  step  is  increased  (see 
Table  1).  Furthermore,  since  is  then  of 

approximately  the  same  value  as  F^"\  the  inner  product 
<F(n+i)*^F(”'''*^>  is  then  also  approaching 
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machine-zero  and  thus,  no  new  direction  vector  can  be 
created  for  the  advancement  of  the  unknown  current  vector. 
The  algorithm  then  becomes  stagnated  and  no  meaningful 
solution  can  be  obtained  no  matter  how  many  more  steps 
are  added. 

A  simple  remedy  which  involves  restarting  the  iterative 
solution  with  a  small  perturbation  of  the  zero  initial 
estimate  —  commonly  adopted  in  most  FFT-based 
implementations —  has  been  suggested  by  Smith,  Peterson 
and  Mittra  [12].  However,  tests  performed  by  the  authors 
usin§  their  ad-^c  procedure  indicated  that  an  enorenous 
solution  may  result  if  insufficient  care  is  taken  in  choosing 
the  amount  of  the  required  perturbation,  although 
convergence  is  usually  alleviated  whenever  the  zero 
starting  vector  is  slightly  perturbed.  The  extent  of 
perturbation,  thus,  plays  an  important  role  in  ensuring  both 
a  convergent  and  correct  solution,  and  since  this  parameter 
is  not  known  in  advance,  gener^  use  of  Smidi  et.  al.'s 
procedure  for  scatterers  other  than  a  simple  conducting 
strip  is  fraught  with  difficulties.  For  higher-dimensiond 
problems,  die  perturbations  in  the  different  coordinate 
components  at  each  discretization  cell  of  the  computational 


domain  would  best  be  related  to  some  known  quantity  of 
the  problem  under  consideration.  At  the  start  of  the  plate 
scattering  simulation,  the  only  known  quantity  is  the 
incident  field  distribution  E'  over  the  plate.  Thus,  it  is 
logical  that  the  initial  estimate  should  be  taken  as  E‘  . 

Fig.  2  shows  the  current  distributions  from  a 
H-polarized  grazing  incidence  achieved  with  this  initial 
estimate  procedure.  Also  included  are  the  CGFFT  solution 
at  two  different  error  tolerances.  It  is  clear  from  this  figure 
that  the  same  behaviour  is  observed  for  this  configuration 
as  for  the  broadside  case  of  Fig.  1,  i.e.  the  Bi-^GFFT 
mediod  requires  a  much  larger  error  tolerances  than  its 
CGFFT  counterpart  when  applied  to  the  same  problem. 
The  numerical  results  associated  with  the  CGFFT  method 
at  the  tolerance  sufficiently  required  by  the  Bi-CGFFT  are, 
again,  not  recognizable  —  particularly  in  the 
y  -  component,  where  the  spikes  at  the  back  comers  are 
certainly  non-physical.  The  final  results  for  the 
Bi-CGFFT  at  e  =  0.3  and  the  CGFFT  at 
€  =0.5*  lO"'  are  virtually  the  same. 


C.  Oblique  incidences 

Next,  the  use  of  E‘  as  an  initial  estimate  for  other  angles 
of  incidence  is  investigated.  Fig.  3  shows  the  convergence 
rates  of  Ae  Bi-CGFFT  method  when  subject  to  two 
arbitrary  oblique  angles  of  incidence.  The  erratic 
behaviours  of  Ae  two  convergence  curves  are  the 
consequences  of  the  fact  that  successive  residual  norms  in 
a  Biconjugate  Gra^ent  formulation  need  not  satisfy  the 

inequality  ||  ||  <  ||  H  at  all  iteration  steps  as 

would  be  expected  in  a  CGFFT  formulation. 


Iteration  number 


(a)  Oi  =  30®  (E-polarized) 


Fig.  3:  Comparison  of  convergence  rates  between  the  two 
methods  of  estimating  the  initial  unknown  from  a 
Bi-CGFFT  solution. 

It  is  worth  emphasizing,  however,  that,  bamiig  the  case  of 
non-convergence  for  certain  grazing  incidences,  the 
Aeoretical  finite-termination  property  of  a  Biconjugate 
Gradient  approach  is  still  obtainable  as  with  a  Conjugate 
Gradient  formulation,  albeit  with  infinite-precision 
arithmetics  [9].  However,  as  Fig.  3  (a)  clearly  shows,  using 
E‘  as  an  initial  estimate  can  substantially  reduced  the 
computational  cost  required  to  achieve  convergence  as 
compared  witii  simply  using  zero.  The  increase  in  iteration 

count  by  using  =  0  in  this  case  is  a  staggering  1680%. 


A  similar  increase  rate  has  also  been  recorded  for  the 
near-grazing  angle  of  89®,  where  although  convergence  is 
achieved,  the  huge  computational  cost  makes  it  necessary 
to  use  E‘  as  an  initial  estimate  instead  of  zero.  On  the  other 
hand,  the  detrimental  effect  of  using  =  E‘  at 
di  =  45®  in  Fig.  3  (b)  is  minimal. 

The  convergence  behaviours  of  the  two  angles  in  Fig.  3 
are  typical  of  a  Bi-CGFFT  solution  when  used  with  two 
different  initial  estimate  procedures  for  other  angles  of 
incidence  in  boA  polarizations;  either  the  convergence 

rates  are  considerably  accelerated  when  =  0  is 

replaced  by  =  E‘,  or  only  neglible  decrease  in 
convergence  rate  is  observed  [14].  In  most  cases,  the 
significant  speedup  ratios  indicate  that  the  incident  electric 
field  E‘  should  be  used  as  an  initial  estimate  for  the 
unknown  current  distribution  in  a  Bi-CGFFT  method. 

D.  The  effect  of  using  the  incident  electric 
field  as  an  initial  estimate  in  a  CGFFT 
formulation 

Using  E'  as  an  initial  estimate  in  a  CGFFT  formulation 
is,  however,  not  beneficial  to  its  resulting  convergence  rate 
as  compart  to  the  usual  procedure  of  simply  using 

=  0.  Fig.  4  shows  the  convergence  rates  associated 
with  three  arbitrary  angles  of  incidence  in  a  CGFFT 
solution.  The  differences  in  convergence  rates  associated 
with  the  two  starting  procedures  appear  to  be  smaller  as  the 
angle  of  incidence  increases.  At  grazing  incidence,  where 

0i  =  90®,  the  discrepancy  is  minimal,  whereas  at 

broadside  incidence  (0,-  =  0®),  it  is  best  to  use  zero  as  a 
starting  estimate.  The  study  in  this  section  shows  clearly 
that  no  advantages  can  be  gained  by  using  E‘  as  an  initial 
estimate  for  the  unknown  current  distribution  in  a  CGFFT 
application  and  zero  distribution  should  be  adopted  as  a 
starting  estimate  whenever  a  CGFFT  approach  is  used  to 
solve  the  plate  scattering  problem. 


Fig.  4:  Comparison  of  convergence  rates  of  the  two  methods 
of  estimating  the  initial  unknown  with  a  CGFFT 
approach  on  a  lAo  *  lAo  plate. 
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I.  CONCLUSIONS 

In  this  paper,  the  t^licability  of  the  Biconjugate 
Gradient  FFT  applied  to  a  pulse-basis  formulation  for  the 
conducting  plate  problem  has  been  demonstrated.  The 
computational  costassociated  with  the  Bi-CGFFT  is  much 
less  than  a  conventional  CGFFT  ^proach  due  to  the 
reduction  of  a  matrix-vector  product  ^r  iteraticMi  step  and 
a  much  larger  error  tolerance  to  achieve  a  similar 
well-behaved  result  as  collared  to  that  obtained  by  its 
CGFFT  counterpart.  The  difncul^  with  non-convergence 
of  the  Bi-CGFPl  at  grazing  incidence  was  allevia^  by 
using  the  incident  electric  field  E‘  as  an  initial  estimate. 
Furthermore,  it  was  also  demonstrated  that  using  this 
starting  guess  procedure  also  accelerates  the  convergence 
rates  associated  with  other  angles  of  incidence. 

The  use  of  E‘  as  an  initial  estimate  is  recently  shown  to 
ovCTCome  the  non-convergence  problem  associated  with 
the  application  of  the  Bi-CGFFT  method  to  2-  and  3-D 
dielectric  bodies  [15].  Although  the  cause  of 
non-convergence  is  different  from  diat  associated  with  a 

thin  conducting  plate,  the  fact  that  using  =  E‘  can 
induce  convergence  is  a  particularly  pleasing  aspect. 
Although  no  justification  is  available  to  explam  for  the 
improvement  in  convergence  rates,  diis  non  ad-hoc  initial 
starting  procedure  appears  to  be  mandatory  in  the 
successful  implementation  of  a  Biconjugate  Gradient  FFT 
method. 

The  Biconjugate  Gradient  FFT  has  been  increasingly 
used  as  an  efficient  algorithm  to  solve  for  the 
discrete-convolutional  system  which  arises  either  ftom  a 
FFT-based  formulation  similar  to  that  considered  in  this 
paper  [16],  or  as  part  of  a  hybrid  method  [17].  The 
mcorporation  of  the  starting  guess  procedure  proposed  by 
the  authors  in  this  paper  may  well  help  to  dleviate  any 
non-convergence  problem  associated  with  the  application 
of  the  Bi-CGFFT  to  otho*  problems  of  interest  and  allay  the 
doubt  concerning  its  consistent  performance  for 
computational  electromagnetics. 
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THE  APPLIED  COMPUTATIONAL  ELECTROMAGNETICS  SOCIETY 
ANNOUNCES  A  SPECIAL  ISSUE  OF  THE  ACES  JOURNAL  ON: 

ADVANCES  IN  THE  APPLICATION  OF  THE 
METHOD  OF  MOMENTS 

TO  ELECTROMAGNETIC  RADIATION  AND  SCATTERING  PROBLEMS 

The  Applied  Computational  Electromagnetics  Society  is  pleased  to  announce  the  publication  of  a  1995 
Special  Issue  of  the  ACES  Journal  on  the  use  of  the  Method  of  Moments  in  the  evaluation  of  electromagnetic 
radiation  and  scattering  problems.  The  objectives  of  this  special  issue  are:  ( 1)  to  provide  the  computational 
electromagnetics  community  with  an  assessment  of  the  current  capabilities  and  uses  of  the  Method  of 
Moments  for  electromagnetics  problems  from  the  low-frequency  to  the  high-frequency  regimes  and  (2)  to 
provide  information  on  recent  advances  that  may  extend  range  of  applicability  and  usefulness  of  the 
Method  of  Moments.  Prospective  authors  are  encouraged  to  submit  papers  of  archival  value  that  address 
these  objectives  and  other  suggested  topics  listed  below. 

SUGGESTED  TOPICS 

•  Modeling  Guidelines  for  Complex  Geometries 

•  Accuracy  Assessment  and  Improvement 

•  Special  Formulations:  Low  Frequency/High  Frequency 

•  Hybrid  Method  of  Moment  Approaches 

•  New  Integral  Equation  Formulations 

•  Peirallelization  of  Moment  Method  Codes 

•  Novel  Equivalence  Principle  Applications  for  the  Method  of  Moments 

•  Fast  Matrix  Computation/Solution  Techniques 

•  Large-Scale  Problems 

DEADLINE  FOR  PAPERS  IS  MARCH  31,  1995 


Send  papers  and  inquiries  to: 


A.W.  Glisson  and  A.A.Kishk 
Special  Guest  Editors 
Department  of  Electrical  Engineering 
University  of  Mississippi 
University,  MS  38677 

Tel:  (601)  232-5353 
FAX:  (601)232-7231 
E-mail:  eeallen@vm.cc.olemiss.edu 
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The  Applied  Computational  Electromagnetics  Society 
Announces  a  Special  Issue  of  the  ACES  Journal  on: 

Applied  Mathematics:  meeting  the  challenges  presented  by 
Computational  Electromagnetics 

The  objectives  of  this  special  issue  are  a)  to  illuminate  some  of  the  current  mathematical 
techniques  in  computational  electromagnetics,  by  a  series  of  review  or  survey  articles,  and 
b)  to  initiate  and  encourage  interaction  between  the  applied  mathematics  community  on 
the  one  hand,  and  the  electrical  engineers  and  physicists  on  the  other.  Papers  submitted 
must  address  mathematical  problems  arising  in  computational  electromagnetics,  and  the 
conclusions  must  have,  moreover,  some  practical  value.  Contact  the  Guest  Editors. 


Suggested  Topics: 

*  Integral  equations  and  integrodifferential  equations 
-k  Eigenfunction  expansions,  both  interior  and  exterior 

*  Selfadjoint,  as  well  as  non-selfadjoint,  operator  approximation 

*  Singularity  expansion  method,  scattering  poles,  natural  modes 

*  Diffraction  and  asymptotics,  application  of  special  functions 

*  Variational  principles,  Galerkin  and  related  methods 

*  Finite  element  methods,  finite  difference  methods 

*  The  radiation  boundary  problem 

*  Solution  of  large  scale  linear  systems 

k  Eigenvalue  estimation,  especially  for  non-selfadjoint  problems 
k  Optimization,  conjugate  and  bi-conjugate  gradient  methods,  GMRES 

*  Numerical  evaluation  of  integrals  with  oscillatory  or  singular  integrands 


Wherever  possible,  attention  should  be  given  to  error  estimates.  This  is  not  just  a  difficult 
mathematical  issue^,  it  is  absolutely  vital  for  all  engineering  considerations  and  it  is  still 
largely  unresolved  in  computational  electromagnetics.  Modern  numerical  analysis  seems  not 
to  have  lived  up  to  its  basic  premise,  because  we  do  not  yet  possess  useful  error  estimates. 


The  deadline  for  papers  is  June  30,  1995. 
Mail  one  hard  copy  to  each  of  the  Guest  Editors: 


Eugene  Tomer 

Applied  Mathematics  and  Computing 

150  Hernandez  Avenue 

San  Francisco,  CA  94127,  USA 

Tel:  (415)  665-9555  Fax:  (415)  731-3551 

e-mail:  etomer@netcom.com 


Andrew  F.  Peterson 

School  of  Electrical  Engineering 

Georgia  Institute  of  Technology 

Atlanta,  GA  30332,  USA 

Tel:  (404)  853-9831  Fax:  (404)  853-9171 

e-mail:  apl6@prism.gatech.edu 


^  See,  e.g.,  S.G.  Mikhlin,  Error  Analysis  in  Numerical  Processes,  Wiley,  1991. 
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PUBLICATION  CRITERIA 

Each  paper  is  required  to  manifest  some  relation  to  applied 
computational  electromagnetics.  Papers  may  address 
general  issues  in  applied  computational  electromagnetics, 
or  they  may  focus  on  speciBc  applications,  techniques, 
codes,  or  computational  issues.  While  the  following  list  is 
not  exhaustive,  each  paper  will  generally  relate  to  at  least  one 
of  these  areas: 

1.  Code  validation.  This  is  done  using  internal  checks  or 
experimental,  analytical  or  other  computational  data. 
Measured  data  of  potential  utility  to  code  validation  efforts 
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identification  of  numerical  accuracy  or  other  limitations, 
solution  convergence,  numerical  and  physical  modeling  error, 
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address  issues  such  as  ease-of-use,  set-up  time,  run  time, 
special  outputs,  or  other  special  features. 

3.  Computational  studies  of  basic  physics.  This  involves 
using  a  code,  algorithm,  or  computational  technique  to 
simulate  reality  in  such  a  way  that  better  or  new  physical 
insight  or  understanding  is  achieved. 

4.  New  computational  techniques,  or  new  applications  for 
existing  computational  techniques  or  codes. 

5.  "Tricks  of  the  trade"  in  selecting  and  applying  codes  and 
techniques. 

6.  New  codes,  algorithms,  code  enhancement,  and  code 
fixes.  This  category  is  self-explanatory  but  includes 
significant  changes  to  existing  codes,  such  as  applicability 
extensions,  algorithm  optimization,  problem  correction, 
limitation  removal,  or  other  performance  improvement. 
NOTE:  CODE  (OR  ALGORITHM)  CAPABILITY 
DESCRIPTIONS  ARE  NOT  ACCEPTABLE,  UNLESS 
THEY  CONTAIN  SUFFICIENT  TECHNICAL  MATERIAL 
TO  JUSTIFY  CONSIDERATION. 

7.  Code  input/output  issues.  This  normally  involves 
innovations  in  input  (such  as  input  geometry  standardization, 
automatic  mesh  generation,  or  computer-aided  design)  or  in 
output  (whether  it  be  tabular,  graphical,  statistical, 
Fourier-transformed,  or  otherwise  signal-processed).  Material 
dealing  with  input^output  database  management,  output 
interpretation,  or  other  input/output  issues  will  also  be 
considered  for  publication. 


8.  Computer  hardware  issues.  This  is  the  category  for 
analysis  of  hardware  capabilities  and  limitations  in  meeting 
various  types  of  electromagnetics  computational  require¬ 
ments.  Vector  and  parallel  computational  techniques  and 
implementation  are  of  particular  interest. 

Applications  of  interest  include,  but  are  not  limited  to, 
antennas  (and  their  electromagnetic  environments),  networks, 
static  fields,  radar  cross  section,  shielding,  radiation  hazards, 
biological  effects,  electromagnetic  pulse  (EMP),  electromag¬ 
netic  interference  (EMI),  electromagnetic  compatibility 
(EMC),  power  transmission,  charge  transport,  dielectric  and 
magnetic  materials,  microwave  components,  MMIC  technol¬ 
ogy,  remote  sensing  and  geophysics,  communications 
systems,  fiber  optics,  plasmas,  particle  accelerators,  genera¬ 
tors  and  motors,  electromagnetic  wave  propagation, 
non-destructive  evaluation,  eddy  currents,  and  inverse 
scattering. 

Techniques  of  interest  include  frequency-domain  and 
time-domain  techniques,  integral  equation  and  differential 
equation  techniques,  diffraction  theories,  physical  optics, 
moment  methods,  finite  differences  and  finite  element 
techniques,  modal  expansions,  perturbation  methods,  and 
hybrid  methods.  This  list  is  not  exhaustive. 

A  unique  feature  of  the  Journal  is  the  publication  of  unsuc¬ 
cessful  efforts  in  applied  computational  electromagnetics. 
Publication  of  such  material  provides  a  means  to  discuss 
problem  areas  in  electromagnetic  modeling.  Material 
representing  an  unsuccessful  application  or  negative  results  in 
computational  electromagnetics  will  be  considered  for 
publication  only  if  a  reasonable  expectation  of  success  (and 
a  reasonable  effort)  are  reflected.  Moreover,  such  material 
must  represent  a  problem  area  of  potential  interest  to  the 
ACES  membership. 

EDITORIAL  REVIEW 

In  order  to  ensure  an  appropriate  level  of  quality  control, 
papers  are  refereed.  They  are  reviewed  both  for  technical 
correctness  and  for  adherence  to  the  listed  guidelines  regard¬ 
ing  information  content.  Authors  should  submit  the  initial 
manuscript  in  draft  form  so  that  any  suggested  changes  can  be 
made  before  the  photo-ready  copy  is  prepared. 
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usable  by  technical  abstracting  and  indexing  services. 
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ence,  we  recommend  that  references  be  given  by  author(s) 
name  and  year  in  the  body  of  the  paper  (with  alphabetical 
listing  of  all  references  at  the  end  of  the  paper).  Titles  of 
Journals,  monographs,  and  similar  publications  should  be  in 
boldface  or  italic  font  or  should  be  underlined.  Titles  of 
papers  or  articles  should  be  in  quotation  marks. 

5.  Internal  consistency  shall  also  be  maintained  for  other 
elements  of  style,  such  as  equation  numbering.  As  a  guideline 
for  authors  who  have  no  other  preference,  we  suggest  that 
equation  numbers  be  placed  in  parentheses  at  the  right 
column  margin. 

6.  The  intent  and  meaning  of  all  text  must  be  clear.  For 
authors  who  are  NOT  masters  of  the  English  language,  the 
ACES  Editorial  Staff  will  provide  assistance  with  grammar 
(subject  to  clarity  of  intent  and  meaning). 

7.  Unused  space  should  be  minimized.  Sections  and  subsec¬ 
tions  should  not  normally  begin  on  a  new  page. 
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mall  address  of  at  least  one  of  the  authors.  Only 
camera-ready  original  copies  are  accepted  for  publication, 
although  authors  may  submit  other  copies  for  publication 


review.  The  term  "camera-ready"  means  that  the  material  is 
neat,  legible,  and  reproducible.  The  preferred  font  style  is 
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While  ACES  reserves  the  right  to  re-type  any  submitted 
material,  this  is  not  generally  done. 
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